分块矩阵的逆

最编程 2024-02-13 20:58:40

...

分块矩阵怎么求逆？ - 知乎四分块矩阵求逆是有通式的（参见分块矩阵求逆公式）下面是一些分块矩阵求逆公式：另可参考Block matrix o…https://www.zhihu.com/question/47760591Matrix inversion identities - 知乎Recently I reviewed some tricks in obtaining the inversion of the matrices. An efficient technique to obtain the inversion of a matrix makes use of Schur's Complement. In China, this trick is al…https://zhuanlan.zhihu.com/p/157436207

已经验证了上面两个公式是正确的。

import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.metrics.pairwise import rbf_kernel
from time import time
from numpy.linalg import inv

X, y = datasets.load_iris(return_X_y=True)

A = rbf_kernel(X=X[0:5],gamma=0.1) + 0.1 * np.eye(5)
B = rbf_kernel(X=X[0:5],Y=[X[5]],gamma=0.1)
C = B.T
D = rbf_kernel(X=[X[5]],Y=[X[5]],gamma=0.1) + 0.1
# print(A.shape)
# print(B.shape)
# print(C.shape)
# print(D.shape)

K = rbf_kernel(X=X[0:6],gamma=0.1) + 0.1 * np.eye(6)
print("直接")
print(np.linalg.inv(K))
# S_time = time()
# for i in range(10000):
#     AAi = np.linalg.inv(K)
# print("Time_1::",time()-S_time)

Ai = np.linalg.inv(A)
S_time = time()
average_time = []
for i in range(10000):
    n = A.shape[0]
    a1 = np.concatenate((np.eye(n),-Ai@B),axis=1)
    b1 = np.concatenate((np.zeros((1,n)),np.eye(1)),axis=1)
    M1 = np.concatenate((a1,b1),axis=0)
    # print("M1:",M1.shape)

    a2 = np.concatenate((Ai,np.zeros((n,1))),axis=1)
    b2 = np.concatenate((np.zeros((1,n)),np.linalg.inv(D-C@Ai@B)),axis=1)
    M2 = np.concatenate((a2,b2),axis=0)
    # print("M2:",M2.shape)

    a3 = np.concatenate((np.eye(n),np.zeros((n,1))),axis=1)
    b3 = np.concatenate((-C@Ai,np.eye(1)),axis=1)
    M3 = np.concatenate((a3,b3),axis=0)
    # print("M2:",M3.shape)

    s_time = time()
    AAi = M1@M2@M3
    e_time = time()
    cost_time = e_time - s_time
    average_time.append(cost_time)
    # print(M1@M2@M3)
# print("Time_2::",time()-S_time)
print("Time_2::",np.mean(average_time))

DCAB = inv(D-C@Ai@B)
# for i in range(10000):
a = Ai+ Ai@B@DCAB@C@Ai
b = -Ai@B@DCAB
c = -DCAB@C@Ai
d = DCAB
M1 = np.concatenate((a,b),axis=1)
M2 = np.concatenate((c,d),axis=1)
M = np.concatenate((M1,M2),axis=0)
print("M:")
print(M)

import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.metrics.pairwise import rbf_kernel
from time import time
from numpy.linalg import inv

X, y = datasets.load_iris(return_X_y=True)

A = rbf_kernel(X=X[0:100],gamma=0.1) + 0.1 * np.eye(100)
B = rbf_kernel(X=X[0:100],Y=[X[100]],gamma=0.1)
C = B.T
D = rbf_kernel(X=[X[100]],Y=[X[100]],gamma=0.1) + 0.1


K = rbf_kernel(X=X[0:101],gamma=0.1) + 0.1 * np.eye(101)
# print("直接")
# print(np.linalg.inv(K))
S_time = time()
for i in range(1000):
    AAi = np.linalg.inv(K)
print("Time_1::",time()-S_time)

Ai = np.linalg.inv(A)
# S_time = time()
# # average_time = []
# for i in range(1000):
#     n = A.shape[0]
#     a1 = np.concatenate((np.eye(n),-Ai@B),axis=1)
#     b1 = np.concatenate((np.zeros((1,n)),np.eye(1)),axis=1)
#     M1 = np.concatenate((a1,b1),axis=0)
#     # print("M1:",M1.shape)
#
#     a2 = np.concatenate((Ai,np.zeros((n,1))),axis=1)
#     b2 = np.concatenate((np.zeros((1,n)),np.linalg.inv(D-C@Ai@B)),axis=1)
#     M2 = np.concatenate((a2,b2),axis=0)
#     # print("M2:",M2.shape)
#
#     a3 = np.concatenate((np.eye(n),np.zeros((n,1))),axis=1)
#     b3 = np.concatenate((-C@Ai,np.eye(1)),axis=1)
#     M3 = np.concatenate((a3,b3),axis=0)
#     # print("M2:",M3.shape)
#
#     s_time = time()
#     AAi = M1@M2@M3
#     e_time = time()
#     cost_time = e_time - s_time
#     # average_time.append(cost_time)
#     # print(M1@M2@M3)
# print("Time_2::",time()-S_time)
# # print("Time_2::",np.mean(average_time))

DCAB = inv(D-C@Ai@B)
S_time = time()
for i in range(1000):
    a = Ai+ Ai@B@DCAB@C@Ai
    b = -Ai@B@DCAB
    c = -DCAB@C@Ai
    d = DCAB
    M = np.concatenate((np.concatenate((a,b),axis=1),np.concatenate((c,d),axis=1)),axis=0)
print("Time_3::",time()-S_time)

S_time = time()
for i in range(1000):
    a = Ai+ Ai@B@DCAB@C@Ai
    b = -Ai@B@DCAB
    c = -DCAB@C@Ai
    d = DCAB
    M1 = np.concatenate((a,b),axis=1)
    M2 = np.concatenate((c,d),axis=1)
    M = np.concatenate((M1,M2),axis=0)
print("Time_4::",time()-S_time)

使用分块矩阵计算缺失节约了大量的计算代价（前提A的逆已知）。

import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.metrics.pairwise import rbf_kernel
from time import time
from numpy.linalg import inv

X, y = datasets.load_iris(return_X_y=True)

A = rbf_kernel(X=X[0:5],gamma=0.1) + 0.1 * np.eye(5)
B = rbf_kernel(X=X[0:5],Y=[X[5]],gamma=0.1)
C = B.T
D = rbf_kernel(X=[X[5]],Y=[X[5]],gamma=0.1) + 0.1


K = rbf_kernel(X=X[0:6],gamma=0.1) + 0.1 * np.eye(6)
print("直接")
print(np.linalg.inv(K))
S_time = time()
for i in range(10000):
    AAi = np.linalg.inv(K)
print("Time_1::",time()-S_time)

n = A.shape[0]
Ai = np.linalg.inv(A)
# S_time = time()
# # average_time = []
# for i in range(1000):
#     n = A.shape[0]
#     a1 = np.concatenate((np.eye(n),-Ai@B),axis=1)
#     b1 = np.concatenate((np.zeros((1,n)),np.eye(1)),axis=1)
#     M1 = np.concatenate((a1,b1),axis=0)
#     # print("M1:",M1.shape)
#
#     a2 = np.concatenate((Ai,np.zeros((n,1))),axis=1)
#     b2 = np.concatenate((np.zeros((1,n)),np.linalg.inv(D-C@Ai@B)),axis=1)
#     M2 = np.concatenate((a2,b2),axis=0)
#     # print("M2:",M2.shape)
#
#     a3 = np.concatenate((np.eye(n),np.zeros((n,1))),axis=1)
#     b3 = np.concatenate((-C@Ai,np.eye(1)),axis=1)
#     M3 = np.concatenate((a3,b3),axis=0)
#     # print("M2:",M3.shape)
#
#     s_time = time()
#     AAi = M1@M2@M3
#     e_time = time()
#     cost_time = e_time - s_time
#     # average_time.append(cost_time)
#     # print(M1@M2@M3)
# print("Time_2::",time()-S_time)
# # print("Time_2::",np.mean(average_time))

DCAB = inv(D-C@Ai@B)
S_time = time()
for i in range(10000):
    a = Ai+ Ai@B@DCAB@C@Ai
    b = -Ai@B@DCAB
    c = -DCAB@C@Ai
    d = DCAB
    M = np.concatenate((np.concatenate((a,b),axis=1),np.concatenate((c,d),axis=1)),axis=0)
print("Time_3::",time()-S_time)

S_time = time()
for i in range(10000):
    a = Ai+ Ai@B@DCAB@C@Ai
    b = -Ai@B@DCAB
    c = -DCAB@C@Ai
    d = DCAB
    M1 = np.concatenate((a,b),axis=1)
    M2 = np.concatenate((c,d),axis=1)
    M = np.concatenate((M1,M2),axis=0)
print("Time_4::",time()-S_time)

S_time = time()
for i in range(10000):
    M = np.zeros((n+1,n+1))
    M[:n,:n] = Ai+ Ai@B@DCAB@C@Ai
    # print(M[:n,n].shape)
    M[:n,n] = -Ai@B@DCAB.reshape(1)
    M[n,:n] = -DCAB@C@Ai
    M[n,n] = DCAB
    # a = Ai+ Ai@B@DCAB@C@Ai
    # b = -Ai@B@DCAB
    # c = -DCAB@C@Ai
    # d = DCAB
    # M1 = np.concatenate((a,b),axis=1)
    # M2 = np.concatenate((c,d),axis=1)
    # M = np.concatenate((M1,M2),axis=0)
print("Time_5::",time()-S_time)

M = np.zeros((n+1,n+1))
M[:n,:n] = Ai+ Ai@B@DCAB@C@Ai
M[:n,n] = -Ai@B@DCAB.reshape(1)
M[n,:n] = -DCAB@C@Ai
M[n,n] = DCAB
print(M)

时间对比（样本数加到100）：

最后一种写法最节约时间。

import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.metrics.pairwise import rbf_kernel
from time import time
from numpy.linalg import inv

X, y = datasets.load_iris(return_X_y=True)

A = rbf_kernel(X=X[0:100],gamma=0.1) + 0.1 * np.eye(100)
B = rbf_kernel(X=X[0:100],Y=[X[100]],gamma=0.1)
C = B.T
D = rbf_kernel(X=[X[100]],Y=[X[100]],gamma=0.1) + 0.1


K = rbf_kernel(X=X[0:101],gamma=0.1) + 0.1 * np.eye(101)
# print("直接")
# print(np.linalg.inv(K))
S_time = time()
for i in range(10000):
    AAi = np.linalg.inv(K)
print("Time_1::",time()-S_time)

n = A.shape[0]
Ai = np.linalg.inv(A)
# S_time = time()
# # average_time = []
# for i in range(1000):
#     n = A.shape[0]
#     a1 = np.concatenate((np.eye(n),-Ai@B),axis=1)
#     b1 = np.concatenate((np.zeros((1,n)),np.eye(1)),axis=1)
#     M1 = np.concatenate((a1,b1),axis=0)
#     # print("M1:",M1.shape)
#
#     a2 = np.concatenate((Ai,np.zeros((n,1))),axis=1)
#     b2 = np.concatenate((np.zeros((1,n)),np.linalg.inv(D-C@Ai@B)),axis=1)
#     M2 = np.concatenate((a2,b2),axis=0)
#     # print("M2:",M2.shape)
#
#     a3 = np.concatenate((np.eye(n),np.zeros((n,1))),axis=1)
#     b3 = np.concatenate((-C@Ai,np.eye(1)),axis=1)
#     M3 = np.concatenate((a3,b3),axis=0)
#     # print("M2:",M3.shape)
#
#     s_time = time()
#     AAi = M1@M2@M3
#     e_time = time()
#     cost_time = e_time - s_time
#     # average_time.append(cost_time)
#     # print(M1@M2@M3)
# print("Time_2::",time()-S_time)
# # print("Time_2::",np.mean(average_time))

DCAB = inv(D-C@Ai@B)
S_time = time()
for i in range(10000):
    a = Ai+ Ai@B@DCAB@C@Ai
    b = -Ai@B@DCAB
    c = -DCAB@C@Ai
    d = DCAB
    M = np.concatenate((np.concatenate((a,b),axis=1),np.concatenate((c,d),axis=1)),axis=0)
print("Time_3::",time()-S_time)

S_time = time()
for i in range(10000):
    a = Ai+ Ai@B@DCAB@C@Ai
    b = -Ai@B@DCAB
    c = -DCAB@C@Ai
    d = DCAB
    M1 = np.concatenate((a,b),axis=1)
    M2 = np.concatenate((c,d),axis=1)
    M = np.concatenate((M1,M2),axis=0)
print("Time_4::",time()-S_time)

S_time = time()
for i in range(10000):
    M = np.zeros((n+1,n+1))
    M[:n,:n] = Ai+ Ai@B@DCAB@C@Ai
    # print(M[:n,n].shape)
    M[:n,n] = -Ai@B@DCAB.reshape(1)
    M[n,:n] = -DCAB@C@Ai
    M[n,n] = DCAB
    # a = Ai+ Ai@B@DCAB@C@Ai
    # b = -Ai@B@DCAB
    # c = -DCAB@C@Ai
    # d = DCAB
    # M1 = np.concatenate((a,b),axis=1)
    # M2 = np.concatenate((c,d),axis=1)
    # M = np.concatenate((M1,M2),axis=0)
print("Time_5::",time()-S_time)

上一篇：块矩阵逆的推导过程

下一篇：理解矩阵求逆的引理和分块求逆的推导过程

分块矩阵的逆

已经验证了上面两个公式是正确的。

使用分块矩阵计算缺失节约了大量的计算代价（前提A的逆已知）。

二叉树, 329. 矩阵中的最长递增路径, 备忘递归表填充, [73] [算法分析与设计] 516.

相关运算、卷积运算和托普利兹矩阵之间的关系

基于 Numpy 的 Python 数组和矩阵_python

教程 | 入门：深度学习矩阵操作的概念和代码实现

三维空间中的刚体运动 - (1) 驰放坐标和旋转矩阵

带有平移缩放旋转功能的 opengl 矩阵变换

[华为 OD] 给定一个 N*M 矩阵，求矩阵每列 M 个元素的最大值 100

矩阵终极大 Boss 的作用是什么？

切线空间 (TBN) ---- 谈论图形中的矩阵运算

矩阵范式的特性和证明