Finding the index k of square nilpotent matrix A

I am trying to find the k index of a nilpotent square matrix A. It doesn't seems to do the job precisely tho..

function k = nilIndex(A)    
       [n,n]=size(A);
       l=0;
       k=n;
       B =A^l;
       while k-l>1
           m=round((l+k)/2);
           d=m-l;
           C=A^(d)*B;
           if all(abs(C) == 0, 'all')
               k = m;
           else 
                   l = m;
                   B = A^l;
           end
       end
       if k-l==1
          k=n; 
       end
end

Answer 1

[John D'Errico]

A nilpotent matrix is a square matrix A that has the property that for SOME integer power k, we have A^k is entirely zero. The smallest power to cause it to be entirely zero is called the index. And for a matrix of size nxn, the index k cannot exceed n. As such, it is probably reasonable to just use a simple loop, where each iteration will just mutiply by A. Check for all zeros at each iteration.

Can you get more fancy? Well yes. If the matrix is large, then it is probably efficient to compute some subset of powers of A. Store them away. For a matrix of size nxn, I would guess that an appropriate scheme would have you generate the powers of A from 1 to sqrt(n), saving all of them. If any of those powers were all zero, then you are done. If not, then you can use a variation of a bisection search to find the actual index. (You just need to get tricky. NEVER raise the matrices to a power.)

Is that more efficient than a purely linear search form 1 to n? ONLY if n is quite large would it even matter.

As for using a direct bisection search, that forces you to compute relatively large powers of A at each iteration. And that in itself will not be efficient. For example, consider the random matrix:

n = 2500;
A = rand(n,n); % surely not nilpotent, but who cares as an example?
timeit(@() A*A)

ans = 0.0741

timeit(@() A^1000)

ans = 0.9050

So raising a large matrix to a power is considerably more expensive than a simple matrix multiply. With n = 2500, it looks like you can do roughly 15 simple matrix multiplies for each such power operation. How would you write it as a loop? DON"T USE POWERS!

As a test, we know a tridiagonal matrix with all zeros on the diagonal must be nilpotent.

A = triu(ones(9),1)    

A = 9x9

0 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

nilIndex(A)

ans = 9

This next matrix will surely not be nilpotent.

A = randi(10,[10,10]) - 5

A = 10x10

-4 2 1 3 -2 5 -4 2 -2 0 5 4 0 0 -4 5 5 -2 -1 4 2 -4 1 -3 -2 4 -2 2 0 0 -1 3 -2 2 4 3 1 2 4 3 1 2 -1 4 -3 -3 3 -3 -3 -4 -1 -3 -1 0 -2 -4 3 1 2 -4 1 0 4 4 4 4 -1 2 4 3 5 4 -4 1 -3 5 -1 2 -2 -4 -4 4 3 -4 4 0 -3 2 -2 4 -3 3 -3 -2 4 -3 2 1 -3 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

nilIndex(A)

A is not nilpotent

ans = NaN

function k = nilIndex(A)    
n = size(A,1);
k = NaN;
Apow = A;
for i = 1:n
  if ~any(Apow,'all')
    k = i;
    break
  else
    Apow = Apow*A;
  end
end
if isnan(k)
  disp('A is not nilpotent')
end
end

The point is, your bisection scheme is not nearly as efficient as you think it is. Again, a simple loop is surely better than you think. So, unless this was a homework assignment, where you were instructed to use bisection (sigh) or you have the coding chutzpah to code up that scheme I described above, just use a basic loop. An interesting question is what would be the optimal size of the initial linear sample to employ.

Answer 2

[Ayush Anand]

Hi,

You can do away with the binary search entirely and do a simple linear search if n is not large:

function k = nilIndex(A)
n = size(A, 1); % Assuming A is square
k = 1; % Start checking from A^1
while k <= n
    if all(all(A^k == 0)) 
        return; % Return current k if A^k is zero
    end
    k = k + 1;
end
k = NaN; % Return NaN if no such k is found within n iterations 
end

However, if you want to stick to the binary search approach, you can get rid of the

if k-l==1
    k=n
end

part as in a scenario where k-l is 1 at the end of the loop, k does not need to be reset to n. An example would be for a matrix like A = [[3,3,3] ; [4,4,4] ; [-7,-7,-7]], at the end of the loop the value of l is 1 and the value of k is 2. In this case k-l = 1 at the end of the loop but k =2 is the correct answer and does not need to be reset to 3.

Hope this helps!

Answer 3

[Bruno Luong]

I modify your original code and now it seems working:

function k = nilIndex(A)
[n,n]=size(A);
l=0;
k=n;
B =A^l;
while k-l>1
    m=round((l+k)/2);
    d=m-l;
    C=A^(d)*B;
    if all(C == 0, 'all') % no need abs
        k = m;
    else
        l = m;
        B = C; % aleady computed
    end
end
% why destroy the binary search?
       %if k-l==1
       %   y search.k=n; 
       %end
end

Answer 4

[Bruno Luong]

This modified version of binary search only use matrix multiplication and keep track of A^2^p, p = 0,1,2 ...The number of matrix products is 2*ceil(log2(k)) It seems the fastest

% Generate random binary nilpotent matrix
A = diag(rand(1,1000)>0.01, 1);
p = randperm(length(A));
A = A(p,p);

tic
k = nilpotentdegree(A);
toc

Elapsed time is 0.217295 seconds.

if isempty(k)
    fprintf('A is not nilpotent\n')
else
    % This test must return 1
    test_k = all(A^k == 0,'all') && any(A^(k-1),'all'),
    fprintf('A nilpotent with degree = %d\n', k)
end

test_k = logical
   1

A nilpotent with degree = 333

function k = nilpotentdegree(A)
% Check if A is nilpotent, return the degree
% k is empty if A is not impotent
n = size(A,1);
if size(A,2) ~= n
    error('A lusr be square')
end

q = nextpow2(n);
A2p = zeros([n n q+1]);
% Compue A^(2^p) for p = 0,1,2, ..., q
% stored in A2p(:,:,l) with l := p+1;
% A^2^q is 0 if A is nilpotent
AP = A;
tp = 1;
for p = 0:q
    if nnz(AP) && tp<n
        A2p(:,:,p+1) = AP;    
        AP = AP*AP;
        tp = tp*2;
    else
        break
    end
end

if nnz(AP)
    k = [];
else
    Ak = eye(n);
    k = 1;
    for l = p:-1:1
        tp = tp/2;
        Bk = Ak * A2p(:,:,l);
        if nnz(Bk)
            k = k + tp;
            Ak = Bk;
        end
    end
end
end