1

I am currently working on a program that aims to take some steady-state probability distribution (specified by particular free parameters of my choosing) and propagates it to some other steady-state (also specified by the same set of modified parameters). I can do this by calculating a tridiagonal matrix (not symmetric, unfortunately) and finding its eigenvalues/eigenvectors, such that:

$$p^{(2)}=\hat{K}^{(2)}(t,t_0)p^{(1)}=V^{(2)}Y^{(2)}U^{(2)}p^{(1)}$$

where the 2's denote the second probability distribution. I let $\hat{K}^{(2)}(t,t_0)$ denote the propagator of the second probability distribution. V is a matrix of the normalized right eigenvectors, U is a matrix of the normalized left eigenvectors, and Y is a diagonal matrix whose entries are not important to our discussion (it is the identity when t=0). The left and right eigenvectors form an orthonormal set and hence I should get that:

$$V^{(2)}U^{(2)}=I$$

I can obtain $U^{(2)}$ from $V^{(2)}$ simply by taking the inverse of $V^{(2)}$, i.e.

$$\big(V^{(2)}\big)^{-1}=U^{(2)}$$

Here is where the issue comes on. $V^{(2)}$ is very close to being a singular matrix. This would be fine if the matrix were not extremely large (300x300-500x500), as MATLAB's Moore-Penrose pseudo-inverse would work beautifully in obtaining the information I need; however, since the matrix is so large, I am not getting results that are necessarily accurate (for example, if I set the parameters that specify the shape of the later probability distribution to a point past a certain threshold, the propagator does not turn into the identity at t=0, as $V^{(2)}U^{(2)} \neq I$ for large portions of the matrix). I believe this error stems from the Moore-Penrose pseudo-inverse algorithm, as it is unable to build a proper inverse of $V^{(2)}$. I am hence curious if there is a way around my predicament - perhaps a better algorithm for calculating pseudo-inverses and/or inverses for matrices that are extremely close to being singular in general.

user107053
  • 33
  • You have to improve the condition number by preconditioning. There are lots of preconditioners for linear systems out there. One thing that would help would be symmetrization provided the eigenvalues turn out to be real. That said, it sounds to me like you should be using the SVD instead. – Ian Jul 11 '17 at 16:17
  • @Ian The pseudo-inverse algorithm that MATLAB uses (pinv(B)) uses the Moore-Penrose algorithm, which actually computes the pseudo-inverse using singular value decomposition. – user107053 Jul 11 '17 at 16:34
  • 1
    I meant the SVD of $K$ itself. You probably already have severe errors by the time you have just computed $V$ and $Y$. – Ian Jul 11 '17 at 16:38
  • If $V$ is very close to singular, its inverse will be very sensitive to roundoff error, and no numerical method for inverting $V$ is likely to give you a good result for $U$. It may be better to compute the left eigenvectors separately by finding the left nullspace of $A - \lambda I$ for each eigenvalue $\lambda$. – Robert Israel Jul 11 '17 at 17:54

0 Answers0