Finding the index k of square nilpotent matrix A
Show older comments
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
Accepted Answer
John D'Errico
on 12 Apr 2024
Edited: John D'Errico
on 12 Apr 2024
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)
timeit(@() A^1000)
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)
nilIndex(A)
This next matrix will surely not be nilpotent.
A = randi(10,[10,10]) - 5
nilIndex(A)
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.
More Answers (3)
Ayush Anand
on 12 Apr 2024
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!
Bruno Luong
on 15 Apr 2024
Edited: Bruno Luong
on 15 Apr 2024
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
Bruno Luong
on 15 Apr 2024
Edited: Bruno Luong
on 15 Apr 2024
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
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
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
Categories
Find more on Matrices and Arrays in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
