CODEWITHSUNDEEP
pythonperformancelinear-algebrasparse-matrixmatrix-factorization
aMATHeur
aMATHeur

Reputation: 41

LDL factorization of a sparse matrix in Python

Let's say I have a sparse symmetric real indefinite matrix A of size N=150,000. This matrix is generated at each iteration of an algorithm with the same structure (positions of non-zero elements) but different values. I need to solve at each iteration several systems of the form Ax=b and to compute the matrix inertia (numbers of positive, negative and zero eigenvalues). From available literature, using an LDL.T factorization of the matrix A seems a promising way.

A is a KKT matrix, of the following form:

|--------+-----|
|        |     |
|   H    | G.T |
|        |     |
|--------+-----|
|   G    |  0  |
|--------+-----|

With H mostly diagonal, and G is rectangular, made of overlapping rectangular blocks along the diagonal (obtained by discretization of a state equations system).

I am using python. I have tried using the numpy package, the cvxopt package for dense and sparse implementations, but the factorization is way too slow for my goal.

Has anybody a python package / method to recommend for such a (medium) size problem? For example, a good implementation of a multifrontal direct solver/factorizer? Is there any way to exploit to the fullest the fact that the structure of the matrix is the same from one iteration to the next? Any thoughts outside of the box?

Upvotes: 3

Views: 587

Answers (1)

Andrino
Andrino

Reputation: 36

Others might have a slighly different problem so first some context. Specifically if A would be symetric positive definite one can simply do a cholesky decomposition and then split off the diagonal. For this numerous software support exists. For quasi-definite matrices this is not possible hence the questioner asks about the LDL decomposition, a generallization of Cholesky.

The QDLDL software package provides an LDL implementation in C. It is part of the well known OSQP quadratic convex optimization solver. https://github.com/osqp/qdldl. There is also a python interface: https://github.com/osqp/qdldl-python that is installable with PIP pip install qdldl

Solving the linear system Ax=b is achievable by:

import qdldl
F = qdldl.Solver(A)
x = F.solve(b)

Butefully QDLDL also allows to only update the specific values in A when the sparsity pattern remains. In these cases the symbolic factorization stage has to be done only once. The computational cost is greatly reduced.

F.update(A_new)

From a research perspective I would appreciate if You are able to provide context about the application and context of the linear system solution.

Upvotes: 0

Related Questions