The LAPACK routines xGETRI compute the inverse of a matrix $A = PLU$ in its LU decomposed form by first computing $U^{-1}$, and then solving the system: $$ (A^{-1} P) L = U^{-1} $$ My question is: wouldn't it be more numerically stable to compute $L^{-1}$ since $L$ is a unit lower triangular matrix, and thus has determinant 1. Then you could solve: $$ U (A^{-1} P) = L^{-1} $$ for the matrix $A^{-1} P$. The matrix $U$ is not unit triangular, so inverting it could lead to problems if it is close to singular. It seems the alternative scheme of inverting $L$ would be more numerically stable, so I am curious if anyone know why LAPACK chose to invert $U$ instead?


I think the reasoning is as follows:

When one stores the LU decomposition, the diagonal (factorized) of the matrix is stored in $U$, that's why it is not unit-triangular. This is a choice, and one could have opted for storing it in $L$.

Now, the fact that the matrix is not-invertible, does not mean that it does not have LU factorization. Take a look at getrf description for the info output parameter:

info=i, $u_{ii}$ is 0. The factorization has been completed, but $U$ is exactly singular. Division by $0$ will occur if you use the factor $U$ for solving a system of linear equations.

So, there, technically, is no error in getrf call resultant in a non-zero positive info parameter. However, the usefulness of the result is questionable, since both the solution of the system of linear equations, as well as the inversion via getri will result in an error.

So, during the inversion via getri, the invertibility is determined by $U$. Hense, the choice.

I do not see, why would you avoid problems when the matrix is close to singular by flipping what is inverted, you probably would just postpone them until the multiplication by $U$ with very small diagonal entries (in case of close-to-singularness) in case of $$ U (A^{-1} P) = L^{-1} $$

compared to solving against a RHS with very large entries in $U^{-1}$ in case of $$ (A^{-1} P) L = U^{-1} $$

Moreover, pivoting is also involved. This is straightforward (assuming $A$ and $U$ are invertible) $$ A=PLU \Rightarrow I=A^{-1}PLU\Rightarrow U^{-1}=A^{-1}PL $$ I am not sure how to get to $L^{-1}$ version (assuming $A$ and $L$ are invertible) without the change of how pivoting $P$ is applied: $$ \require{cancel} \begin{aligned} A=PLU \Rightarrow I=PLUA^{-1} &\Rightarrow \cancel{?}\\ &\Rightarrow P=PLUA^{-1}P\Rightarrow I=LUA^{-1}P\Rightarrow L^{-1}=UA^{-1}P \end{aligned} $$

Thus, pivoting $P$ has to change to: $$ \require{cancel}\cancel{A=LU\tilde{P}\Rightarrow I=LU\tilde{P}A^{-1}\Rightarrow L^{-1}=U\tilde{P}A^{-1}} $$

which would imply a lot of other issues without obvious benefits.

| cite | improve this answer | |
  • $\begingroup$ Interesting and you are probably correct. I just want to point out you don't need to change the pivoting for the $L^{-1}$ version: $I = PLUA^{-1} \Rightarrow P = PLUA^{-1}P \Rightarrow I = LUA^{-1}P \Rightarrow L^{-1} = UA^{-1} P$. $\endgroup$ – vibe Jun 10 '19 at 21:15
  • $\begingroup$ @vibe does not it has to assume that $P$ is invertible (step 2 in your derivation in the comments)? $\endgroup$ – Anton Menshov Jun 10 '19 at 21:17
  • $\begingroup$ @vibe and now I have a question to myself if $P$ is always invertible (which it should be....) :) $\endgroup$ – Anton Menshov Jun 10 '19 at 21:18
  • 2
    $\begingroup$ The determinant of a permutation matrix is always 1 or -1, so they are invertible. In fact they form a group under matrix multiplication en.m.wikipedia.org/wiki/Permutation_matrix $\endgroup$ – vibe Jun 10 '19 at 21:56

Let $A = PLU$ be an LU decomposition, where $A$ is $n\times n$. Note that $U = D(I+N)$ where $D$ is diagonal, $I$ is the identity matrix, and $N$ is strictly upper-triangular. Now note that $(I+N)^{-1} = I + \sum_{k = 1}^{n-1} (-1)^k N^k$, and so $U^{-1} = \left(I + \sum_{k = 1}^{n-1} (-1)^k N^k \right) D^{-1}$, which can be computed pretty efficiently it you look at it hard enough. You can invert $L$ the same way, the only difference is you don't scale the columns of the "factorization" by $D^{-1}$. So it seems to me that the two procedures are pretty much equivalent, it just depends on if you want to scale by $D^{-1}$ at the start when "inverting" your triangular matrix, or at the end when solving the system.

| cite | improve this answer | |
  • $\begingroup$ I found a reference: Du Croz and Higham, Stability of methods for matrix inversion, IMA J. Num. Anal., 1992. They analyze several methods of inverting a matrix, including $A^{-1} = U^{-1} L^{-1}$, $A^{-1} L = U^{-1}$, and a couple others, and find they are all roughly equivalent in terms of numerical stability. Interestingly, LAPACK chose the $A^{-1} L = U^{-1}$ method, which requires extra workspace, while the $U^{-1} L^{-1}$ method requires no extra workspace. $\endgroup$ – vibe Jun 18 '19 at 17:29

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.