我需要反转一个大而密集的矩阵,我希望使用Scipy的
这是我的问题:为了编写计算矩阵向量积
MATLAB有类似的设置,其中我们编写
我的问题是,我无法弄清楚如何让
谢谢!
gmres
来完成。幸运的是,密集矩阵A遵循一种模式,我不需要在内存中存储矩阵。 LinearOperator
类允许我们构建一个对象,它可以作为GMRES的矩阵,并直接计算矩阵向量积A*v
。也就是说,我们编写一个函数mv(v)
,它以向量v
作为输入,并返回mv(v) = A*v
。然后,我们可以使用LinearOperator
类创建A_LinOp = LinearOperator(shape = shape, matvec = mv)
。我们可以将线性算子放入Scipy的gmres
命令中,以评估矩阵向量乘积,而无需完全加载A
到内存中。
LinearOperator
的文档在这里:LinearOperator
Documentation。这是我的问题:为了编写计算矩阵向量积
mv(v) = A*v
的例程,我需要另一个输入向量C
。 A
中的条目形式为A[i,j] = f(C[i] - C[j])
。因此,我真正想要的是,mv
有两个输入,一个是固定向量输入C
,另一个是可变输入v
,我们希望计算A*v
。MATLAB有类似的设置,其中我们编写
x = gmres(@(v) mv(v,C),b)
,其中b
是问题Ax = b
的右手边,mv
是以可变输入v
作为变量的函数,我们希望计算A*v
,而C
是已知的固定向量,我们需要用它来组装A
。我的问题是,我无法弄清楚如何让
LinearOperator
类接受两个输入,一个可变和一个“固定”,就像在MATLAB中一样。是否有一种方法在SciPy中执行类似的操作?或者,如果有人知道更好的方法来反转一个大而密集的矩阵(50000, 50000)
,其中条目遵循一种模式,我将非常感激任何建议。谢谢!
编辑:我应该明确这个信息。该矩阵实际上是(分块形式)[A C; C^T 0]
,其中A
是N x N
(N
很大),C
是N x 3
,0
是3 x 3
,而C^T
是C
的转置。这个数组C
与上面提到的那个数组是相同的。A
的条目遵循模式A[i,j] = f(C[i]-C[j])
。
我写了mv(v,C)
来逐行构造A*v[i]
,其中i=0,N
,通过计算求和f(C[i]-C[j)*v[j]
(实际上我做的是numpy.dot(FC,v)
,其中FC[j] = f(C[i]-C[j])
,这个方法效果不错)。然后,在结束时完成对C^T
行的计算。我希望最终用multiprocessing
调用来并行化for循环,但这是考虑未来的事情。我还将研究使用Cython加速计算。
functools.partial
)替代它:make_mv = lambda c: lambda v: sum(f(c[i] - c[j])...
(但为了可读性最好使用def
),而且这个函数可能应该用Cython编写,否则对于这么大的矩阵来说速度会太慢。 - jorgecamv2
函数将使用当前作用域中的任何C
,而不是创建函数时所在作用域中的C
。如果你想让mv2
引用特定的C
,你绝对需要上述两种方法之一,然后使用mv = make_mv(C)
。 - jorgeca