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!

  • $\begingroup$ 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$. $\endgroup$ – Wolfgang Bangerth May 30 at 17:11
  • $\begingroup$ 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. $\endgroup$ – LAD May 30 at 17:14
  • $\begingroup$ @WolfgangBangerth I have edited my post to make it clearer. $\endgroup$ – LAD May 30 at 17:22
  • 2
    $\begingroup$ 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 $\endgroup$ – vibe May 31 at 0:54
  • $\begingroup$ This also looks like you might want to read through the literature on sparse direct solvers. $\endgroup$ – 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.

| cite | improve this answer | |
  • $\begingroup$ Thank you so much for your detailed answer! I will read it carefully. $\endgroup$ – LAD Jun 1 at 13:12
  • 1
    $\begingroup$ I should have mentioned, a C implementation of recursive LU can be found here: git.savannah.gnu.org/cgit/gsl.git/tree/linalg/lu.c -- see the routine LU_decomp_L3. Of course, you don't really need a fully recursive algorithm, you are just doing the first step of the recursion, taking advantage of the special matrix structure, and then calling a conventional LU algorithm $\endgroup$ – vibe Jun 1 at 16:50
  • $\begingroup$ I would like to thank you again for your dedicated help. For a student not specializing in coding and algorithm. Adjusting available algorithm is really exhausting :((( $\endgroup$ – LAD Jun 1 at 16:58
  • $\begingroup$ 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. $\endgroup$ – vibe Jun 1 at 17:01
  • $\begingroup$ My professor allows me to choose Python or C, C++, or C#. I will change to C, C++, or C# if it makes my life easier. $\endgroup$ – LAD Jun 1 at 17:02

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.