I need to compute the Drazin inverse $A^D$ of a singular M-matrix $A$, i.e., a matrix in the form $A=\lambda I -P$, where $P$ has nonnegative entries and $\lambda$ is the spectral radius (Perron value) of $P$. I already know the right and left eigenvectors $v$ and $u^T$ of $P$ with eigenvalue $\lambda$ (that is, the vectors in the left and right kernel of $A$).

1) Is there a Matlab subroutine around for computing Drazin inverses? I can't seem to find any, so I had to create my own (which is probably very inefficient)

2) Is there a way to exploit the knowledge of the two nullspaces (and the fact that they are 1-dimensional) to speed up this computation?

  • $\begingroup$ Have you seen dx.doi.org/10.1137/S0895479891228279 already? $\endgroup$ – J. M. is not a mathematician Oct 2 '10 at 1:32
  • $\begingroup$ @J. M. I don't have access to the full-text from my institution, but from the abstract it does not seem related to what I'm doing. I already know the Perron right and left eigenvectors, I need to compute the group/Drazin inverse $\endgroup$ – Federico Poloni Oct 2 '10 at 12:34
  • $\begingroup$ By the way, an answer to question (ii) for the Moore-Penrose pseudo-inverse instead of the Drazin inverse would be most welcome as well. $\endgroup$ – Federico Poloni Oct 2 '10 at 12:35
  • $\begingroup$ Hey, I don't know if you've already found an answer to this, but it looks to me that the results of dx.doi.org/10.1007/BFb0120751 and dx.doi.org/10.1137/0131057 might be applicable to your problem. $\endgroup$ – J. M. is not a mathematician Dec 14 '11 at 13:28
  • $\begingroup$ The first link is a "DOI not found". I'll check the second, though, thanks! $\endgroup$ – Federico Poloni Dec 14 '11 at 20:11

If you have access to Matlab, and the matrix is not gigantic you can use the following to calculate the Drazin inverse:

  1. Determine the index of $A$, i.e, does $\mathrm{rank}(A^k) = \mathrm{rank}(A^{k+1})$

  2. Solve the linear system $A^{k+1} X = A^k$

  3. Use the fact that $A^D$, the Drazin inverse, is given by

$$A^D = A^k X^{k+1}$$

Jim Shoaf, N. C. Central University

  • $\begingroup$ Thanks! Could you please expand on how to solve that singular linear system? If I am not missing some crucial detail, it could even be $0X=0$, for instance in the case when A=[0 1; 0 0]. $\endgroup$ – Federico Poloni Aug 26 '15 at 22:06
  • $\begingroup$ 1) One way to try to solve A^(k+1)*X = A^k in Matlab is X = A^(k+1) \ A^k , where these powers of A are stored as Ak2 and Ak 2) Another approach to calculating Ad is the Schulzz (Hoteling) iteration: X = .0001*A repeat X = 2*X - XAX (until convergence -- sometimes divergence!!) $\endgroup$ – Jim Shoaf Jul 25 at 4:41

This one might have some pertinent leads:



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.