在Python中实现MINRES

Question

有没有可以处理厄米特矩阵的MINRES伪逆算法的python实现？

我找到了一些源码，但所有这些源码都只能处理实矩阵，而且似乎不容易推广到复杂的情况：

https://searchcode.com/codesearch/view/89958680/

https://github.com/pascanur/theano_optimize

(还有其他几个链接，但我的声誉不允许我发布它们)

Answer 1

一个大小为$n$的厄米特系统

$$\mathbf y= \mathbf H^{-1}\mathbf v$$

可以嵌入到大小为$2n$的真实对称系统中：

\begin{等式} \begin{bmatrix} \Re(\mathbf y)\Im(\mathbf y) \end{bmatrix} = \begin{bmatrix} \Re(\mathbf H)&-\Im(\mathbf H)\Im(\mathbf H)&\Re(\mathbf H) \end{bmatrix}^{-1} \begin{bmatrix} \Re(\mathbf v)\Im(\mathbf v) \end{bmatrix}。\end{公式}

最小残差方法通常用于大型问题，其中构造$H$是不切实际的。在这种情况下，我们可以有一个运算来计算矩阵向量乘积$f：\mathbb C^n \to \mathbb C^n；，，f(\mathbf v) = \mathbf H\mathbf v$这个函数可以被包装为在$\mathbf x \in \mathbb R^{2n}$上操作，方法是将$\mathbf x$转换回复向量，应用$f$，然后将结果嵌入到$\mathbb R^{2n}$中。

以下是python / numpy /scipy中的示例：

from scipy.sparse.linalg import minres, LinearOperator
from pylab import *

# Problem size
N = 100

# error helper
er = lambda t,a,b:print('%s error:'%t,mean(abs(a-b)))

# random Hermitian matrix
Q = randn(N,N) + 1j*randn(N,N)
H = Q@conj(Q.T)

# random complex vector
v = randn(N) + 1j*randn(N)

# ground-truth solution
x0 = inv(H)@v

# Pack/unpack complex vector as stacked real vector
c2r = lambda v:block([real(v),imag(v)])
r2c = lambda v:kron([1,1j],eye(N))@v

# Verify that we can embed C^n in R^(2N)
Hr = real(H)
Hi = imag(H)
Hs = block([[Hr,-Hi],[Hi,Hr]])
vs = c2r(v)
xs = inv(Hs)@vs
x1 = r2c(xs)
er('Embed',x0,x1)

# Verify that minres works as expected in R-embed
x2 = r2c(minres(Hs,vs,tol=1e-12)[0])
er('Minres 1',x0,x2)

# Demonstrate using operators
Av = lambda u:c2r( H @ r2c(u) )
A  = LinearOperator((N*2,)*2,Av,Av)

# Minres, converting input/output to/from complex/real
x3 = r2c(minres(Hs,vs,tol=1e-12)[0])
er('Minres 2',x0,x3)

>>> Embed error: 5.317184726020268e-12
>>> Minres 1 error: 6.641342200989796e-11
>>> Minres 2 error: 6.641342200989796e-11