用奇异块反转块对角矩阵

Question

我有一个块对角矩阵，其中一些对角块为零。矩阵是一个运算符的表示，我想要这个运算符的逆（各个块的条件足够低，我不用担心）

import numpy as np

np.random.seed(42)

N = 10

def composite_operator():
    matrix = np.zeros((3*N, 3*N))
    for n in range(1,N):
        matrix[3*n:3*n+3, 3*n:3*n+3] = np.random.rand(3,3)
    return matrix

上面的代码块应该有助于重新创建与我现在所拥有的类似的东西，除了 N~1000 并且在我的例子中，有多个零对角线块沿对角线分布在已知位置。

与该运算符相反，我希望零值块仍然为零。有没有一种方法可以真正有效且快速地做到这一点？目前我正在创建一个新的

composite_operator_inv

并用除了零块之外的倒置块一一填充它，我真的很希望有一种更快的方法来做到这一点。

Answer 1

您可以创建一个 3 x 3 块的矩阵，并按以下方式对其求逆。我添加了执行时间测量来显示该方法的效率：

import numpy as np
import scipy.linalg as la
import time

np.random.seed(42)
N = 1000  

def composite_operator(N):
    start_time = time.time() 
    blocks = []
    for n in range(N):
        if np.random.rand() > 0.1:  
            blocks.append(np.random.rand(3,3))
        else:
            blocks.append(np.zeros((3,3)))  
    matrix = la.block_diag(*blocks)
    end_time = time.time()  
    print("Time to create the matrix: {:.6f} seconds".format(end_time - start_time))
    return matrix

matrix = composite_operator(N)

def invert_composite_operator(matrix):
    start_time = time.time()  
    N = matrix.shape[0] // 3
    inv_blocks = []
    for n in range(N):
        block = matrix[3*n:3*n+3, 3*n:3*n+3]
        if np.all(block == 0):
            inv_blocks.append(np.zeros((3,3)))
        else:
            inv_blocks.append(la.inv(block))
    inverse_matrix = la.block_diag(*inv_blocks)
    end_time = time.time()  
    print("Time to invert the matrix: {:.6f} seconds".format(end_time - start_time))
    return inverse_matrix

inverse_matrix = invert_composite_operator(matrix)

与

Time to create the matrix: 0.025877 seconds
Time to invert the matrix: 0.057063 seconds

有了

N = 10000

，你还有合理的时间：

Time to create the matrix: 0.334444 seconds
Time to invert the matrix: 0.424029 seconds