如何使用LU分解求解指定随机生成的50×50矩阵的逆矩阵

Got it, let's walk through this problem step by step—first getting the matrix generation right, then computing its inverse using LU decomposition.

First, a quick note on your reference code: np.random.randint(2) generates a matrix of 0s and 1s (integers). Depending on whether you need floating-point values in [0,1] or binary 0/1 integers, use one of these snippets (both use your specified seed 1007092020):

• For [0,1) floating-point matrix:

import numpy as np
np.random.seed(1007092020)
A = np.random.rand(50, 50)  # Values range from 0 (inclusive) to 1 (exclusive)

• For binary 0/1 integer matrix (matching your reference code's logic):

import numpy as np
np.random.seed(1007092020)
A = np.random.randint(2, size=(50, 50))  # Only 0 and 1 as elements

LU factorization breaks down matrix A into a lower triangular matrix L and upper triangular matrix U (plus a permutation matrix P for numerical stability), so P @ A = L @ U. To find A⁻¹, we use the fact that A @ A⁻¹ = I (the identity matrix). This translates to solving two triangular systems:

1. First solve L @ Y = P.T @ I using forward substitution
2. Then solve U @ A⁻¹ = Y using backward substitution

Option 1: Use SciPy's Built-in Tools (Recommended for Real-World Use) #

SciPy has optimized functions for LU decomposition and triangular system solving—this is the most reliable approach:

from scipy.linalg import lu, solve_triangular

# Perform LU decomposition with permutation matrix
P, L, U = lu(A)

# Step 1: Solve L @ Y = P.T (since P @ A = L @ U, we adjust the identity matrix equation)
Y = solve_triangular(L, P.T, lower=True)

# Step 2: Solve U @ A_inv = Y to get the inverse
A_inv = solve_triangular(U, Y, lower=False)

# Verify the result: A @ A_inv should be nearly the identity matrix
print(np.allclose(A @ A_inv, np.eye(50)))  # Should print True if correct

Option 2: Manual Implementation (For Learning Purposes) #

If you want to understand the underlying math, here's a simplified manual implementation (note: this assumes A is non-singular and doesn't require row swaps—use the SciPy version for robustness):

def lu_decomposition(mat):
    n = mat.shape[0]
    L = np.eye(n)
    U = mat.copy()
    
    for col in range(n - 1):
        for row in range(col + 1, n):
            factor = U[row, col] / U[col, col]
            L[row, col] = factor
            U[row, col:] -= factor * U[col, col:]
    return L, U

def forward_substitution(L, b):
    n = L.shape[0]
    y = np.zeros_like(b)
    for i in range(n):
        y[i] = (b[i] - np.dot(L[i, :i], y[:i])) / L[i, i]
    return y

def backward_substitution(U, y):
    n = U.shape[0]
    x = np.zeros_like(y)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
    return x

# Decompose the matrix
L, U = lu_decomposition(A)

# Solve for each column of the inverse (since each column of I gives a column of A⁻¹)
A_inv = np.zeros_like(A)
for i in range(50):
    b = np.eye(50)[:, i]
    y = forward_substitution(L, b)
    x = backward_substitution(U, y)
    A_inv[:, i] = x

# Verify correctness
print(np.allclose(A @ A_inv, np.eye(50)))

Key Notes #

• Invertibility Check: Not all matrices are invertible! For your 0/1 matrix, first check if its rank is 50 with np.linalg.matrix_rank(A)—if not, it's singular and has no inverse.
• Numerical Stability: The manual implementation skips row permutations, which can lead to division by zero or large errors. Always use the permutation-aware SciPy method for real applications.

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