病态对称矩阵的逆

计算科学 矩阵 精确
2021-12-05 14:08:02

我有一个矩阵K,其维度为,其中每个元素都使用以下等式计算:(n,n)

Ki,j=exp(αti2γ(titj)2αtj2)

...其中内等距分布数字序列。这个矩阵是对称的,所以我希望它的逆矩阵也是对称的。t(n)[3,3]

由于元素趋于零的速度有多快,因此反转该矩阵很困难,但是如果向对角线添加一个小的正数,则基数 R 和 numpy 设法反转矩阵。

问题是计算的这个逆不是对称的。


我认为这可能是由于精度问题。此外,从 ^{-1}中减去的转置会产生一些相当大的值(这是有道理的 - 如果您在 K 中有非常小的值您会期望在),但这会导致问题。有没有办法通过使用一些特殊的库来计算正确的逆矩阵(即对称的并且在再次反转时只要 K 保持对称,我也不介意对K1K1KK1KK


编辑:Sympy 支持任意精度的矩阵求逆,但下面答案中的 vec 技巧和关于 matlab 的 inv 函数的评论非常有趣。

3个回答

,您也不需要逆如果您对这个表达式感兴趣,我将解释如何将其转换为矩阵方程,然后更有效地求解:K1hhTK1K1

让我们将定义为:X

X=K1hhTK1K1

您的目标是通过假设都是已知XKhhT

KXK=hhTK

我定义:

B=hhTK

通过使用运算符,该方程将转换为标准线性方程:vec

(KTK)vec(X)=vec(B)

定义:KTK=A

最后:

Avec(X)=vec(B)

有许多有效的线性求解器,正如您在评论中提到的的维数不是那么高,它应该是相当有效的。K

更新:

此 Python 代码将建议的方法与直接计算的方法进行比较:X

import numpy as np

d = 10

K = np.zeros((d,d))

alpha = 1
gamma = 10

t_vec = np.linspace(-3,3,d)

for i in range(d):
        for j in range(d):
                K[i][j] = np.exp(-alpha*(t_vec[i]**2)-alpha*(t_vec[j]**2)-gamma*((t_vec[i]-t_vec[j])**2))

A = np.kron(K.T,K)

h = np.random.rand(d,d)

B = np.matmul(h,h.T) - K

vecB = B.flatten(order='F')

vecX = np.linalg.solve(A,vecB)

X = vecX.reshape((d,d))

print X

Kinv = np.linalg.inv(K)

X_direct = np.matmul(np.matmul(Kinv,np.matmul(h,h.T)),Kinv) - Kinv

print X_direct

print np.abs(X-X_direct)

输出是:

X


[[1.59909954e+16 1.82624715e+12 4.24902134e+10 4.53728217e+08
  2.05584992e+08 2.44120203e+08 5.98264287e+08 3.43515584e+10
  3.12705250e+12 1.31655477e+16]
 [1.82624715e+12 5.09588697e+09 4.78071076e+06 2.16534890e+05
  9.59540287e+04 2.92820805e+04 2.33641869e+05 1.22037867e+07
  2.91201539e+09 3.74298362e+12]
 [4.24902134e+10 4.78071076e+06 1.65187655e+05 2.15188820e+03
  8.22687923e+02 8.22528897e+02 2.14605722e+03 1.18784215e+05
  1.27880035e+07 3.05868993e+10]
 [4.53728217e+08 2.16534890e+05 2.15188820e+03 3.48340497e+01
  1.33228000e+01 1.04249016e+01 3.45272712e+01 2.08247428e+03
  2.62962350e+05 3.37679580e+08]
 [2.05584992e+08 9.59540287e+04 8.22687923e+02 1.33228000e+01
  4.33655487e+00 4.31336472e+00 1.42180232e+01 7.12777365e+02
  9.01300211e+04 2.01206353e+08]
 [2.44120203e+08 2.92820805e+04 8.22528897e+02 1.04249016e+01
  4.31336472e+00 4.81918386e+00 1.18107131e+01 7.69453357e+02
  7.36892858e+04 2.10911516e+08]
 [5.98264287e+08 2.33641869e+05 2.14605722e+03 3.45272712e+01
  1.42180232e+01 1.18107131e+01 4.00277702e+01 1.87646704e+03
  1.95000703e+05 5.95359066e+08]
 [3.43515584e+10 1.22037867e+07 1.18784215e+05 2.08247428e+03
  7.12777365e+02 7.69453357e+02 1.87646704e+03 1.41229229e+05
  1.68594796e+07 2.65035020e+10]
 [3.12705250e+12 2.91201539e+09 1.27880035e+07 2.62962350e+05
  9.01300211e+04 7.36892858e+04 1.95000703e+05 1.68594796e+07
  3.77226961e+09 2.08638514e+12]
 [1.31655477e+16 3.74298362e+12 3.05868993e+10 3.37679580e+08
  2.01206353e+08 2.10911516e+08 5.95359066e+08 2.65035020e+10
  2.08638514e+12 1.53843211e+16]]

X_direct

[[1.59909954e+16 1.82624715e+12 4.24902134e+10 4.53728217e+08
  2.05584992e+08 2.44120203e+08 5.98264287e+08 3.43515584e+10
  3.12705250e+12 1.31655477e+16]
 [1.82624715e+12 5.09588697e+09 4.78071076e+06 2.16534890e+05
  9.59540287e+04 2.92820805e+04 2.33641869e+05 1.22037867e+07
  2.91201539e+09 3.74298362e+12]
 [4.24902134e+10 4.78071076e+06 1.65187655e+05 2.15188820e+03
  8.22687923e+02 8.22528897e+02 2.14605722e+03 1.18784215e+05
  1.27880035e+07 3.05868993e+10]
 [4.53728217e+08 2.16534890e+05 2.15188820e+03 3.48340497e+01
  1.33228000e+01 1.04249016e+01 3.45272712e+01 2.08247428e+03
  2.62962350e+05 3.37679580e+08]
 [2.05584992e+08 9.59540287e+04 8.22687923e+02 1.33228000e+01
  4.33655487e+00 4.31336472e+00 1.42180232e+01 7.12777365e+02
  9.01300211e+04 2.01206353e+08]
 [2.44120203e+08 2.92820805e+04 8.22528897e+02 1.04249016e+01
  4.31336472e+00 4.81918386e+00 1.18107131e+01 7.69453357e+02
  7.36892858e+04 2.10911516e+08]
 [5.98264287e+08 2.33641869e+05 2.14605722e+03 3.45272712e+01
  1.42180232e+01 1.18107131e+01 4.00277702e+01 1.87646704e+03
  1.95000703e+05 5.95359066e+08]
 [3.43515584e+10 1.22037867e+07 1.18784215e+05 2.08247428e+03
  7.12777365e+02 7.69453357e+02 1.87646704e+03 1.41229229e+05
  1.68594796e+07 2.65035020e+10]
 [3.12705250e+12 2.91201539e+09 1.27880035e+07 2.62962350e+05
  9.01300211e+04 7.36892858e+04 1.95000703e+05 1.68594796e+07
  3.77226961e+09 2.08638514e+12]
 [1.31655477e+16 3.74298362e+12 3.05868993e+10 3.37679580e+08
  2.01206353e+08 2.10911516e+08 5.95359066e+08 2.65035020e+10
  2.08638514e+12 1.53843211e+16]]

np.abs(X - X_direct)

[[4.00000000e+00 4.88281250e-04 1.52587891e-05 0.00000000e+00
  2.98023224e-08 0.00000000e+00 0.00000000e+00 1.14440918e-05
  4.88281250e-04 6.00000000e+00]
 [1.46484375e-03 3.81469727e-06 1.02445483e-08 8.73114914e-11
  2.91038305e-11 2.54658516e-11 2.91038305e-11 5.58793545e-09
  9.53674316e-07 9.76562500e-04]
 [0.00000000e+00 3.72529030e-09 5.82076609e-11 0.00000000e+00
  3.41060513e-13 3.41060513e-13 9.09494702e-13 5.82076609e-11
  1.86264515e-09 3.81469727e-06]
 [5.96046448e-08 1.45519152e-10 4.54747351e-13 2.13162821e-14
  7.10542736e-15 1.77635684e-15 1.42108547e-14 0.00000000e+00
  1.16415322e-10 5.96046448e-08]
 [2.98023224e-08 2.91038305e-11 1.13686838e-13 7.10542736e-15
  1.77635684e-15 3.55271368e-15 3.55271368e-15 2.27373675e-13
  4.36557457e-11 2.98023224e-08]
 [1.19209290e-07 2.91038305e-11 2.27373675e-13 1.77635684e-15
  4.44089210e-15 1.77635684e-15 1.77635684e-15 1.13686838e-13
  4.36557457e-11 2.98023224e-08]
 [0.00000000e+00 2.91038305e-11 4.54747351e-13 0.00000000e+00
  3.55271368e-15 1.77635684e-15 7.10542736e-15 0.00000000e+00
  5.82076609e-11 2.38418579e-07]
 [7.62939453e-06 5.58793545e-09 1.45519152e-11 9.09494702e-13
  1.13686838e-13 2.27373675e-13 4.54747351e-13 5.82076609e-11
  3.72529030e-09 0.00000000e+00]
 [0.00000000e+00 9.53674316e-07 3.72529030e-09 1.16415322e-10
  1.45519152e-11 4.36557457e-11 8.73114914e-11 3.72529030e-09
  1.90734863e-06 2.44140625e-04]
 [2.00000000e+00 4.88281250e-04 7.62939453e-06 1.19209290e-07
  5.96046448e-08 2.98023224e-08 2.38418579e-07 3.81469727e-06
  1.70898438e-03 2.00000000e+00]]

您可以看到差异非常小,并且表明所提出的方法确实有效。

基于Federico 的建议和想法,更直接的提取公式K1hhTK1K1将是:

X=K1hhTK1K1

KXK=hhTK

Z=XK

求解Z

KZ=hhTK

然后找到XT从:

KTXT=ZT

最后X

X=(XT)T

让我们将矩阵的直接求逆与我最初提出的方法和 Federico 的方法之间的误差定义为:

ε=||XdirectX||F||Xdirect||F

在哪里||||F是弗罗贝尼乌斯范数。

我基于 numpy 编写了这段代码:

import numpy as np
import time

d = 100

K = np.zeros((d,d))

alpha = 1
gamma = 10

t_vec = np.linspace(-3,3,d)

for i in range(d):
        for j in range(d):
                K[i][j] = np.exp(-alpha*(t_vec[i]**2)-alpha*(t_vec[j]**2)-gamma*((t_vec[i]-t_vec[j])**2))

A = np.kron(K.T,K)

h = np.random.rand(d,d)

B = np.matmul(h,h.T) - K

vecB = B.flatten(order='F')

start = time.time()

vecX = np.linalg.solve(A,vecB)

end = time.time()

print "My method time elapsed: " + str(end - start) + " seconds"

X = vecX.reshape((d,d))

Kinv = np.linalg.inv(K)

X_direct = np.matmul(np.matmul(Kinv,np.matmul(h,h.T)),Kinv) - Kinv

epsilon1 = np.linalg.norm(X-X_direct) / np.linalg.norm(X_direct)

start = time.time()

Z = np.linalg.solve(K,B)

X = np.linalg.solve(K.T,Z.T).T

end = time.time()

epsilon2 = np.linalg.norm(X-X_direct) / np.linalg.norm(X_direct)

print "Federico's method time elapsed: " + str(end - start) + " seconds"

print "My method error: " + str(epsilon1)

print "Federico's method error: " + str(epsilon2)

这些是结果:

My method time elapsed: 19.122369051 seconds
Federico's method time elapsed: 0.000936031341553 seconds
My method error: 0.99999999989649
Federico's method error: 0.6635256191409429

您会看到 Federico 的方法快了大约 5 个数量级,而且它的误差大约是我提出的方法的一半。

我的第一次尝试是拉出对角线缩放:

Ki,j=exp(αti2γ(titj)2αtj2)=exp(αti2)exp(γ(titj)2)exp(αtj2),
所以K=DMD, 在哪里D是对角线Dii=exp(αti2))Mij=exp(γ(titj)2). 然后你可以重新排列总和并减少到反转M代替K. 您可以尝试扩大正方形并拉出条款exp(γti2)也是,但我不确定生成的矩阵是否更容易反转,因为按照我的建议做M在对角线上有一个,在外面有快速衰减的元素,看起来像一个不错的结构。

我没有尝试过实际做实验,更重要的是你没有指定的值αγ这是重要的信息,但我想这会带来巨大的改进。

[编辑:我已经快速检查了结果矩阵的条件α=1,γ=10; 它从大约减少。1021到约1018, 数字; 所以看起来情况有所改善,但没有太大改善。(而且我无论如何也不相信这些数字,因为它们接近机器精度的倒数。)不过,更重要的是,如果我没记错的话M是一个已知矩阵,即所谓的高斯 Toeplitz 矩阵它的 Cholesky 因式分解有一个封闭形式的表达式,可以帮助您反转它。]

另一个可能有帮助的重新安排是X=K1hhTK1K1=K1(hhTK)K1,它允许您使用linsolve和 LDL 分解而不是显式inverses(您应该知道,这几乎总是一个坏主意)。

我建议你把这两个技巧结合起来。