# Algorithm to factorize matrix whose many rows are already of upper triangular form?

I have a matrix whose many rows are already in the upper triangular form.

$$\begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} & x_{5} \\ 0 & x_{22} & x_{23} & x_{24} & x_{25} \\ 0 & 0 & x_{23} & x_{34} & x_{35} \\ 0 & 0 & 0 & 0 & x_{45} \\ x_{51} & x_{52} & x_{53} & x_{54} & x_{55} \\ x_{61} & x_{62} & x_{63} & x_{64} & x_{65} \end{bmatrix}$$

Let me take an example to show that the function scipy.linalg.lu from package Scipy does not take advantage of this special structure. Here B is a copy of A with the elements below the main diagonal zeroed.

import numpy as np
import scipy.linalg as la
import time
A = np.random.randint(100, size=(10000, 10000))
B = np.triu(A, 0)

start = time.time()
(P, L, U) = la.lu(A)
end = time.time()
print('Time to decompose A =', end - start)

start = time.time()
(P1, L1, U1) = la.lu(B)
end = time.time()
print('Time to decompose B =', end - start)


The result is

Time to decompose A = 5.622066497802734
Time to decompose B = 5.322663068771362


Because my square matrix is of very large dimension and this procedure is repeated thousands of times. I would like to ask for a method (or references) to make use of this special structure to reduce the computational complexity.

Thank you so much for your help!

• Why would you want to decompose a matrix that is already in triangular form? If $A$ is lower triangular, then the $LU$ decomposition is $A=AI$ where $L=A$, $U=I$. – Wolfgang Bangerth May 30 at 17:11
• Hi @WolfgangBangerth, I take this special example in Python to illustrate that scipy.linalg.lu does not make use of this special structure. In my problem, the matrix is not of upper triangular form, but it contains some rows of this form. I mean how to make use of the special structure of these rows. – LAD May 30 at 17:14
• @WolfgangBangerth I have edited my post to make it clearer. – LAD May 30 at 17:22
• gnu.org/software/gsl/doc/html/… This algorithm uses a specialized QR method to factor this type of matrix, namely a modified Elmroth and Gustafson 2000 method – vibe May 31 at 0:54
• This also looks like you might want to read through the literature on sparse direct solvers. – Wolfgang Bangerth May 31 at 2:31

I believe you can accomplish what you want efficiently using the recursive LU algorithm. In brief, recursive LU on a $$M \times N$$ matrix $$A$$ proceeds by partitioning the matrix into 4 blocks: \begin{align} \pmatrix{A_{11} & A_{12} \\ A_{21} & A_{22}} &= \pmatrix{L_{11} & 0 \\ L_{21} & L_{22}} \pmatrix{U_{11} & U_{12} \\ 0 & U_{22}} \\ &= \pmatrix{L_{11} U_{11} & L_{11} U_{12} \\ L_{21} U_{11} & L_{21} U_{12} + L_{22} U_{22}} \end{align} Then, the following 4 subproblems must be solved:

(1) $$A_{11} = L_{11} U_{11}$$ (recursive LU call)

(2) $$A_{12} = L_{11} U_{12} \rightarrow U_{12} = L_{11}^{-1} A_{12}$$ (TRSM - Level 3 BLAS)

(3) $$A_{21} = L_{21} U_{11} \rightarrow L_{21} = A_{21} U_{11}^{-1}$$ (TRSM - Level 3 BLAS)

(4) $$A_{22} = L_{21} U_{12} + L_{22} U_{22} \rightarrow A_{22} - L_{21} U_{12} = L_{22} U_{22}$$ (GEMM followed by recursive LU call)

Now normally for recursive LU, you pick the partition so that $$A_{11}$$ has $$N/2$$ columns. However, for your specialized matrix, you should pick $$A_{11}$$ to be a square upper triangular matrix. So basically just pick $$A_{11}$$ so that you cut off the trapezoidal part of the matrix, and then pick $$A_{12}$$ to be that extra rectangle completing the trapezoid. If you do this, then you immediately know that $$L_{11} = I$$ and $$U_{11} = A_{11}$$ and $$U_{12} = A_{12}$$. Then all that remains to do is steps 3 and 4:

(3) $$L_{21} = A_{21} U_{11}^{-1}$$ (TRSM)

(4) $$A_{22} - L_{21} U_{12} = L_{22} U_{22}$$ (GEMM and then use any standard LU routine, e.g. LAPACK or scipy or whatever)

If you want to also implement pivoting, then it is possible to modify the above steps to allow pivoting.

• You mentioned you are using python, which I am not too familiar with. But for step 3, it is important that you use the proper TRSM BLAS call to take advantage of the triangular structure of $U_{11}$. Hopefully python will allow you to call the appropriate low-level BLAS routines. – vibe Jun 1 at 17:01