필터 지우기
필터 지우기

Finding the index k of square nilpotent matrix A

조회 수: 29 (최근 30일)
c.m.c.m
c.m.c.m 2024년 4월 12일
편집: Bruno Luong 2024년 4월 15일
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

채택된 답변

John D'Errico
John D'Errico 2024년 4월 12일
편집: John D'Errico 2024년 4월 12일
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.

추가 답변 (3개)

Ayush Anand
Ayush Anand 2024년 4월 12일
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
Bruno Luong 2024년 4월 15일
편집: Bruno Luong 2024년 4월 15일
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
Bruno Luong 2024년 4월 15일
편집: Bruno Luong 2024년 4월 15일
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

카테고리

Help CenterFile Exchange에서 Logical에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by