# Fast inversion of a triangular matrix

I need to inverse a matrix $A$ given its $QR$ decomposition. It's a numerical task.

It is told that the inversion should be "possibly cheap". But it does not look like I can do something more efficient than computing $R^{-1}$ and multiplying (if one needs a matrix, not a product of two matrices) it by $Q^{T}$. However, even if a product is OK (no multiplication needed), I'm still ineffective, because I don't know how to inverse a triangular matrix faster than in cubic time.

Are there any tricky algorithms doing that (I'm not asking about something like fast matrix multiplication, it's a stupid homework) or the task only sounds wise but all I have to do is to invert a triangular matrix as it is taught in linear algebra course or using back substitution?

EDIT: I emphasise that the task is to invert a matrix, not to find a solution of a linear system.

-
The best I can come up with is Divisions: $\small n$, Multiplications: $\small 1/6 n^3 +1/2 n^2 - 2/3 n$ . Anyway. R.P.Brent has done much work in optimizing algorithms for formal powerseries, possibly there is something similar about matrices and invresion. Some of his articles are online; what I've found of him is really sophisitcated - perhaps you find something relevant using google with "brent complexity matrices" ... – Gottfried Helms Dec 17 '11 at 4:46
I've found one article where instead of $\small O(x^3)$ the authors work out a characteristic of $\small O(x^{\log_2(7)})$. I've just glanced over the article so far, possibly it is valid only for a subclass of matrices. I've found this online at jstor in "Triangular Factorization and Inversion by Fast Matrix Multiplication", James R. Bunch and John E. Hopcroft Mathematics of Computation Vol. 28, No. 125 (Jan., 1974) (pp. 231-236) – Gottfried Helms Dec 17 '11 at 23:26
Just found this more basic explanative article at "Ask Dr.Math" mathforum.org/library/drmath/view/51908.html which also links further to some research article. – Gottfried Helms Dec 17 '11 at 23:54
Thank you very much! The article by Bunch & Hopcroft is very valuable. You should have written it as an answer. – savick01 Dec 18 '11 at 23:10
That article @Gottfried linked you to is freely available. – J. M. Dec 19 '11 at 0:02

Any $N \times N$ triangular system can be solved in $\mathcal{O}(N^2)$.

For instance, if it is an upper triangular system, start from the last equation $(N^{th}$ equation) which requires only one division to get the $N^{th}$ unknown. Once you have this go to the previous equation $((N-1)^{th}$ equation) which requires only one multiplication, one subtraction and one division to get the $(N-1)^{th}$ unknown. Go to the $(N-2)^{nd}$ equation which requires $2$ multiplications, $2$ subtractions and $1$ division. In general, the $(N-k)^{th}$ equation requires $k$ multiplications, $k$ subtractions and $1$ division. Hence, the total cost is $$1+2+3+\cdots+(N-1) = \frac{N(N-1)}{2} \text{ multiplications}$$ $$1+2+3+\cdots+(N-1) = \frac{N(N-1)}{2} \text{ subtractions}$$ $$N \text{ divisions}$$

The same idea works for lower triangular systems as well (in which case you start from the first equation and proceed all the way down).

In fact, the idea behind all matrix decomposition algorithms is to make the solving part cheaper so that given a linear system even if the right hand side were to change you can solve it in a relatively inexpensive way once you have the decomposition of the matrix. The two main decomposition algorithms which are used satisfy this requirement.

1. $A=LU$. Factoring $A$ into a lower triangular times an upper triangular. The factorization cost is $\mathcal{O}(N^3)$. But once this is done solving $Ax = b$ requires solving $Ly = b$ and $UX = y$ both costing $\mathcal{O}(N^2)$ since both are triangular systems.
2. $A = QR$. Factoring $A$ into an orthonormal matrix times an upper triangular matrix. The factorization cost is $\mathcal{O}(N^3)$. But once this is done solving $Ax = b$ requires solving $Qy = b$ and $RX = y$. The nice thing about orthonormal matrices is that the inverse is nothing but the transpose. Hence, $y=Q^Tb$ which is nothing but a matrix vector product which costs $\mathcal{O}(N^2)$. Solving $Rx = y$ costs $\mathcal{O}(N^2)$ since it is a triangular system.

Other decomposition algorithms like the SVD where $A = U \Sigma V^T$ where $U$ and $V$ are orthonormal and $\Sigma$ is a diagonal also satisfy the requirement. Once we have $A = U \Sigma V^T$, solving $Ax = b$ is equivalent to solving $Uy = b$, whose solution is given by $y = U^Tb$ and costs $\mathcal{O}(N^2)$, $\Sigma z = y$, which can be easily inverted since $\Sigma$ is just a diagonal matrix and hence costs $\mathcal{O}(N)$, and $V^Tx = z$, whose solution is given by $x = Vz$ and costs $\mathcal{O}(N^2)$.

EDIT In case you want the inverse of the lower triangular operator, you proceed as follows. (In numerical linear algebra it is one of the cardinal sins to find the inverse explicitly. In any application you will never need to find the inverse explicitly.) $$L = \begin{pmatrix}1 & 0 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 & 0 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & 0 & \cdots & 0 \\ l_{41} & l_{42} & l_{43} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & l_{n4} & \cdots & l_{nn} \\ \end{pmatrix}$$ We can write $L$ as $$L = L_1 L_2 L_3 \cdots L_{n-1}$$ where $$L_k = \begin{pmatrix}1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & l_{k+1,k} & 1 & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & l_{k+2,k} & 0 & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & l_{k+3,k} & 0 & \cdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & \vdots & 0 & \cdots & \cdots & \ddots & \vdots\\ 0 & 0 & \cdots & l_{n,k} & 0 & \cdots & \cdots & \cdots & 1 \end{pmatrix}$$ Then $$L^{-1} = L_{n-1}^{-1}L_{n-2}^{-1}L_{n-3}^{-1} \cdots L_{1}^{-1}$$ where $$L_k^{-1} = \begin{pmatrix}1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & -l_{k+1,k} & 1 & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & -l_{k+2,k} & 0 & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & -l_{k+3,k} & 0 & \cdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & \vdots & 0 & \cdots & \cdots & \ddots & \vdots\\ 0 & 0 & \cdots & -l_{n,k} & 0 & \cdots & \cdots & \cdots & 1 \end{pmatrix}$$

-
The name for the process is "Gaussian elimination". That may be some help to readers trying to understand how to actually carry it out. – Keith Irwin Dec 16 '11 at 18:05
Yes, that is what I meant by "using back substitution" - it can be used to solve equation $Rx_i=e_i$ and than say that $R^{-1}=[x_1 \ldots x_n]$. It works in time $O(n^2)+O((n-1)^2) + \ldots + O(1^2) = O(n^3)$ just as I said before. I'm asking if there exists something faster. Again: I want to find an inversion of a matrix, not to solve a triangular system. – savick01 Dec 16 '11 at 18:15
A very nice lecture and I thank you for your effort but that is not what I'm asking about (and in fact I know those things). – savick01 Dec 16 '11 at 18:24
@savick01: Why do you need the inverse? Typically in almost applications you never deal with the inverse. You resort to just solving the linear system. If you want the inverse of a unit lower triangular system, you do it as follows. Write $L = L_1 L_2 L_3 \ldots L_n$ where $L_k$'s are unit lower triangular and the only non-zero entries of $L_k$ aprt from the diagonal are the entires in the $k^{th}$ column which is same as the $k^{th}$ column of $L$. $L_k^{-1}$ is a unit lower triangular matrix with non-zero entries in the $k^{th}$ column where the entries are negative of the entries of $L_k$. – user17762 Dec 16 '11 at 18:27
What is wrong about my question? Somebody TOLD me to find an inverse. I can do it in $O(n^3)$ and know that it is not done frequently in applications and I know about solving linear systems, BUT I don't know if the inversion can be done faster than in $O(n^3)$. So I ask about it. – savick01 Dec 16 '11 at 19:03

I agree with Sivaram's assessment that an actual matrix inversion is almost never needed (except in some applications, like forming the variance-covariance matrix in statistics). That being said, there is an $O(n^3)$ method to invert a triangular matrix in place (but note that it takes less effort than the inversion of a general matrix). Pete Stewart shows the lower-triangular version in his book (the version given there is for lower triangular matrices; I hope the needed modifications for the upper triangular version are transparent to you), and there is a FORTRAN implementation in LAPACK.

As a final note, since you said that this is part of a QR decomposition: if you used Householder matrices for generating the orthogonal factor, you should know that it is usually much better to keep the components of the Householder vectors around than the multiplied-out orthogonal matrix. (The usual storage format for QR decompositions is to use the upper triangle of the original for storing the triangular factor, and the lower triangle (and also possibly an auxiliary array) for the Householder vectors that form the orthogonal factor.) See Golub and Van Loan, for instance.

-
Correct me if I am wrong. What Pete Stewart does is the same as evaluating this $L_n^{-1} L_{n-1}^{-1} \cdots L_1^{-1}$ right which costs $n^2$ since there is a tremendous sparsity which can be exploited? – user17762 Dec 17 '11 at 3:29
Not exactly the same; if you look at the text, what's being done is to partition the lower triangular matrix and repeatedly apply forward elimination on the trailing submatrix. – J. M. Dec 17 '11 at 3:36
I don't get a few things: 1) @SivaramAmbikasaran first said that I shall never multiply $L_i^{-1}$ out and now says that there are so much sparsity that it can be done in $O(n^2)$ (so why not if it is so cheap?). 2) Is there really so much sparsity? As I previously wrote, I can't see so much of it. 3) @ J.M. - Stewart's algorithm looks like a simple back substitution or whatever it is called in English. He computes $x_i$ satisfying $x_i^T R=e_i^T$ one by one. There are two nested loops and multiplied vectors are on average linear in size, so we get $O(n^3)$. – savick01 Dec 18 '11 at 22:52
Well, in presence of the article mentioned by @GottfriedHelms I don't expect that we know an algorithm working in square time, especially such a simple one. – savick01 Dec 18 '11 at 23:13
@savick01: Mine was a question to J.M. I was wondering if in the sparsity of $L_k^{-1}$ helped in reducing the cost to $n^2$ instead of $n^3$. – user17762 Dec 18 '11 at 23:28