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

Question

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!

Answer 1

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.