Python中稀疏实对称矩阵惯性快速计算及Mx=b求解方法问询

Great question—since you need to solve Mx = b before computing the inertia, we can kill two birds with one stone using a sparse symmetric factorization that handles both tasks efficiently, without ever converting your matrix to dense form. Here's a practical, Python-based approach tailored to your 100 ≤ N ≤ 1e5 scale:

Why Sparse LDL^T Decomposition Is Perfect for Your Use Case #

For real symmetric matrices, the inertia (count of positive/negative/zero eigenvalues) directly corresponds to the number of positive, negative, and zero diagonal elements in the matrix's LDL^T factorization (thanks to Sylvester's Inertia Theorem). Better yet, this factorization is exactly what we need to solve linear systems like Mx = b efficiently.

Unlike dense eigenvalue solvers (like numpy.linalg.eigvalsh), sparse LDL^T decomposition runs in O(nnz(M)) time (where nnz is the number of non-zero elements in M), which is orders of magnitude faster for large sparse matrices. It also lets you reuse the factorization for repeated solves or inertia calculations, which fits your need for repeated computations.

Step-by-Step Python Implementation #

We'll use scipy.sparse and its SuperLU backend, which natively supports sparse symmetric matrix factorization and exposes both the linear solve and inertia directly.

1. Setup & Input Prep #

First, ensure your matrix is stored in a sparse format optimized for factorization (CSC format is best for this; convert with sp.csc_matrix(M) if needed):

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import splu

2. Perform Sparse Symmetric Factorization #

Run the factorization with symmetric=True to enable optimized symmetric matrix handling:

# Assume M is a scipy.sparse.csc_matrix (real, symmetric)
# b is a numpy array of shape (N,) or (N, 1)
lu = splu(M, symmetric=True)

3. Solve Mx = b #

Reuse the factorization to solve the linear system (you can call lu.solve(b) repeatedly for different b vectors without re-factorizing):

x = lu.solve(b)

4. Extract the Inertia #

The SuperLU factorization object directly provides the inertia as a tuple (positive_count, negative_count, zero_count):

inertia = lu.inertia
print(f"Inertia: +={inertia[0]}, -={inertia[1]}, 0={inertia[2]}")

Key Optimizations & Notes #

• Reuse the Factorization: If you need to solve multiple Mx = b problems or compute inertia multiple times, keep the lu object—re-factoring is the most expensive step, and reusing it cuts down on repeated work.
• Storage Format: Always use CSC format for factorization/solves (CSR is better for matrix-vector multiplies, but CSC is faster for decomposition).
• Numerical Tolerance: If you need to adjust how zero eigenvalues are detected (e.g., for near-singular matrices), you can manually extract the diagonal elements of the D matrix and apply a custom threshold:

  # Extract D matrix diagonal elements (from LDL^T decomposition)
  d = lu.U.diagonal()  # For symmetric factorization, U corresponds to L^T * D
  # Apply custom tolerance to count zeros
  tol = 1e-10
  positive = np.sum(d > tol)
  negative = np.sum(d < -tol)
  zero = np.sum(np.abs(d) <= tol)
  custom_inertia = (positive, negative, zero)
  

• Handling Singular Matrices: The factorization will still work for singular matrices (with zero diagonal elements in D), and the inertia will correctly count the zero eigenvalues.

内容的提问来源于stack exchange，提问作者aMATHeur