MATLAB矩阵幂算法答案

【问题标题】：MATLAB matrix power algorithmMATLAB矩阵幂算法
【发布时间】：2017-10-05 03:33:01
【问题描述】：

我希望将一种算法从 MATLAB 移植到 Python。所述算法的一个步骤涉及采用A^(-1/2)，其中 A 是一个 9x9 方形复矩阵。据我了解，矩阵的平方根（以及它们的逆矩阵）不是唯一的。

我一直在尝试使用 scipy.linalg.fractional_matrix_power 和使用 A^(-1/2) = exp((-1/2)*log(A)) 与 numpy 内置的 expm 和 logm 函数的近似值。前者特别差，只提供小数点后 3 位的精度，而后者对于左上角的元素来说是正确的，但随着您向下和向右移动，它会变得越来越差。这可能是也可能不是表达式的完全有效的数学解决方案，但对于此应用程序来说还不够。

因此，我希望在 Python 中直接实现 MATLAB 的矩阵幂算法，以便每次都能 100% 确认相同的结果。有没有人对这将如何工作有任何见解或文档？该算法的并行性越高越好，因为最终目标是在 OpenCL 中重写它以实现 GPU 加速。

编辑：按要求提供 MCVE：

[[(0.591557294607941+4.33680868994202e-19j), (-0.219707725574605-0.35810724986609j), (-0.121305654177909+0.244558388829046j), (0.155552026648172-0.0180264818714123j), (-0.0537690384136066-0.0630740244116577j), (-0.0107526931263697+0.0397896274845627j), (0.0182892503609312-0.00653264433724856j), (-0.00710188853532244-0.0050445035279044j), (-2.20414002823034e-05+0.00373184532662288j)], [(-0.219707725574605+0.35810724986609j), (0.312038814492119+2.16840434497101e-19j), (-0.109433401402399-0.174379997015402j), (-0.0503362231078033+0.108510948023091j), (0.0631826956936223-0.00992931123813742j), (-0.0219902325360141-0.0233215237172002j), (-0.00314837555001163+0.0148621558916679j), (0.00630295247506065-0.00266790359447072j), (-0.00249343102520442-0.00156160619280611j)], [(-0.121305654177909-0.244558388829046j), (-0.109433401402399+0.174379997015402j), (0.136649392858215-1.76182853028894e-19j), (-0.0434623984527311-0.0669251299161109j), (-0.0168737559719828+0.0393768358149159j), (0.0211288536117387-0.00417146769324491j), (-0.00734306979471257-0.00712443264825166j), (-0.000742681625102133+0.00455752452374196j), (0.00179068247786595-0.000862706240042082j)], [(0.155552026648172+0.0180264818714123j), (-0.0503362231078033-0.108510948023091j), (-0.0434623984527311+0.0669251299161109j), (0.0467980890488569+5.14996031930615e-19j), (-0.0140208255975664-0.0209483313237692j), (-0.00472995448413803+0.0117916398375124j), (0.00589653974090387-0.00134198920550751j), (-0.00202109265416585-0.00184021636458858j), (-0.000150793859056431+0.00116822322464066j)], [(-0.0537690384136066+0.0630740244116577j), (0.0631826956936223+0.00992931123813742j), (-0.0168737559719828-0.0393768358149159j), (-0.0140208255975664+0.0209483313237692j), (0.0136137125669776-2.03287907341032e-20j), (-0.00387854073283377-0.0056769786724813j), (-0.0011741038702424+0.00306007798625676j), (0.00144000687517355-0.000355251914809693j), (-0.000481433965262789-0.00042129815655098j)], [(-0.0107526931263697-0.0397896274845627j), (-0.0219902325360141+0.0233215237172002j), (0.0211288536117387+0.00417146769324491j), (-0.00472995448413803-0.0117916398375124j), (-0.00387854073283377+0.0056769786724813j), (0.00347771689075251+8.21621958836671e-20j), (-0.000944046302699304-0.00136521328407881j), (-0.00026318475762475+0.000704212317211994j), (0.00031422288569727-8.10033316327328e-05j)], [(0.0182892503609312+0.00653264433724856j), (-0.00314837555001163-0.0148621558916679j), (-0.00734306979471257+0.00712443264825166j), (0.00589653974090387+0.00134198920550751j), (-0.0011741038702424-0.00306007798625676j), (-0.000944046302699304+0.00136521328407881j), (0.000792908166233942-7.41153828847513e-21j), (-0.00020531962049495-0.000294952695922854j), (-5.36226164765808e-05+0.000145645628243286j)], [(-0.00710188853532244+0.00504450352790439j), (0.00630295247506065+0.00266790359447072j), (-0.000742681625102133-0.00455752452374196j), (-0.00202109265416585+0.00184021636458858j), (0.00144000687517355+0.000355251914809693j), (-0.00026318475762475-0.000704212317211994j), (-0.00020531962049495+0.000294952695922854j), (0.000162971629601464-5.39321759384574e-22j), (-4.03304806590714e-05-5.77159110863666e-05j)], [(-2.20414002823034e-05-0.00373184532662288j), (-0.00249343102520442+0.00156160619280611j), (0.00179068247786595+0.000862706240042082j), (-0.000150793859056431-0.00116822322464066j), (-0.000481433965262789+0.00042129815655098j), (0.00031422288569727+8.10033316327328e-05j), (-5.36226164765808e-05-0.000145645628243286j), (-4.03304806590714e-05+5.77159110863666e-05j), (3.04302590501313e-05-4.10281583826302e-22j)]]

【问题讨论】：

“精确到小数点后三位”是什么意思？您的问题是 matlab 和 numpy 结果之间的区别吗？只要两个矩阵都产生B^2==A，它们都是有效的结果。您的应用程序不应该依赖任何唯一性（如果没有唯一性），或者提供额外的约束以使答案唯一。
正如我将运行 scipy.linalg.fractional_matrix_power(A,(-1/2)) 并获得一些类似 (-8.50000000e+01 - 1.65625000e+02j) 的元素，并且没有理由让第一个答案正好是 85。问题在应用程序中，该算法得到迭代，以便稍后的结果收敛。我已经对我的 python 端口进行了超过 1000 次迭代的测试，在某些情况下它们仍然存在 10% 左右的错误。
Here 是一个很好的评论。

标签： python matlab scipy matrix-multiplication

【解决方案1】：

我能想到两种解释，在这两种情况下我都指责用户错误。按时间顺序：

理论 #1（微妙的）

我怀疑您将输入矩阵的打印值从一个代码复制到另一个代码中。 IE。当您切换代码时，您会丢弃双精度，这在逆平方根计算过程中会被放大。

作为证据，我将 MATLAB 的平方根平方根与您在 python 中使用的函数进行了比较。出于尺寸考虑，我将展示一个 3x3 的示例，但是——剧透警告——我对 9x9 随机矩阵做了同样的事情，得到了两个条件号为 11.245754109790719 (MATLAB) 和 11.245754109790818 (numpy) 的结果。这应该可以告诉您结果的相似性，而无需保存和加载两个代码之间的实际矩阵。不过我建议你这样做：关键字是scipy.io.loadmat 和savemat。

我所做的是在 python 中生成随机数据（因为这是我更喜欢的）：

>>> import numpy as np
>>> print((np.random.rand(3,3) + 1j*np.random.rand(3,3)).tolist())
[[(0.8404782758300281+0.29389006737780765j), (0.741574080512219+0.7944606900644321j), (0.12788250870304718+0.37304665786925073j)], [(0.8583402784463595+0.13952117266781894j), (0.2138809231406249+0.6233427148017449j), (0.7276466404131303+0.6480559739625379j)], [(0.1784816129006297+0.72452362541158j), (0.2870462766764591+0.8891190037142521j), (0.0980355896905617+0.03022344706473823j)]]

通过将相同的截断输出复制到两个代码中，我保证输入的对应关系。

MATLAB 中的示例：

>> M = [[(0.8404782758300281+0.29389006737780765j), (0.741574080512219+0.7944606900644321j), (0.12788250870304718+0.37304665786925073j)]; [(0.8583402784463595+0.13952117266781894j), (0.2138809231406249+0.6233427148017449j), (0.7276466404131303+0.6480559739625379j)]; [(0.1784816129006297+0.72452362541158j), (0.2870462766764591+0.8891190037142521j), (0.0980355896905617+0.03022344706473823j)]];
>> A = M^(-0.5);
>> format long
>> disp(A)
  0.922112307438377 + 0.919346397931976i  0.108620882045523 - 0.649850434897895i -0.778737740194425 - 0.320654127149988i
 -0.423384022626231 - 0.842737730824859i  0.592015668030645 + 0.661682656423866i  0.529361991464903 - 0.388343838121371i
 -0.550789874427422 + 0.021129515921025i  0.472026152514446 - 0.502143106675176i  0.942976466768961 + 0.141839849623673i

>> cond(A)

ans =

   3.429368520364765

python中的示例：

>>> M = [[(0.8404782758300281+0.29389006737780765j), (0.741574080512219+0.7944606900644321j), (0.12788250870304718+0.37304665786925073j)], [(0.8583402784463595+0.13952117266781894j), (0.2138809231406249+0.6233427148017449j), (0.7276466404
... 131303+0.6480559739625379j)], [(0.1784816129006297+0.72452362541158j), (0.2870462766764591+0.8891190037142521j), (0.0980355896905617+0.03022344706473823j)]]

>>> A = fractional_matrix_power(M,-0.5)

>>> print(A)
[[ 0.92211231+0.9193464j   0.10862088-0.64985043j -0.77873774-0.32065413j]
 [-0.42338402-0.84273773j  0.59201567+0.66168266j  0.52936199-0.38834384j]
 [-0.55078987+0.02112952j  0.47202615-0.50214311j  0.94297647+0.14183985j]]

>>> np.linalg.cond(A)
3.4293685203647408

我的怀疑是，如果您将scipy.io.loadmat 矩阵放入 python 中，进行计算，scipy.io.savemat 结果并使用 MATLAB 将其加载回，您将看到少于 1e-12 的绝对错误（希望更少）结果之间。

理论 #2（面部一）

我怀疑您使用的是 python 2，而您的 -1/2-powered 除法是一个简单的逆运算：

>>> # python 3 below
>>> # python 3's // is python 2's /, i.e. integer division
>>> 1/2
0.5
>>> 1//2
0
>>> -1/2
-0.5
>>> -1//2
-1

所以如果你使用的是 python 2，那么调用

fractional_matrix_power(M,-1/2)

实际上是M 的倒数。显而易见的解决方案是切换到 python 3。不太明显的解决方案是继续使用 python 2（你不应该这样做，如上例所示），但使用

from __future__ import division

在您的每个源文件之上。这将覆盖简单的/ 除法运算符的行为，使其反映python 3 版本，并且您会少一个头痛。

【讨论】：

首先，我实际上并没有从 MATLAB 导入数据。原始算法以 MATLAB 代码的形式存在，我试图从头开始在 python 中重新创建它。两者都从相同的数据文件中读取。在此之前的所有其他步骤（包括我们过去讨论过的一些复杂的广播操作）都保持完全相同的答案。我已经确认fractional_matrix_power 对我来说在 -0.5 下不起作用，但在 0.5 下效果更好，然后取反。这个结果更接近expm(-0.5*logm(A)) 近似值。然而，它仍有高达 10% 的错误率。
我也在用 3.5.2
@Mike 可能相关的是 numpy/scipy 版本；当我在笔记本电脑前时，我会告诉你我的。同时，请确保差异不是由于计算输入的细微差异造成的，即如果您到目前为止还没有做上述测试，请使用一个肯定相同的矩阵。你的矩阵有多奇异/病态？如果较小的子矩阵重现了您的问题，则应将其添加到您的问题中以进行交叉检查（两个程序）。
我正在运行 numpy 1.11.1 和 scipy 0.17.1 来自 MATLAB 的条件数是 11115584961.2。来自 python 的条件数是 1894696172.73。因此，如果我正确理解了这个指标，我的猜测是非常病态的。
Numpy 1.12.0 和 scipy 0.18.1。这里。但是你的问题是：这些条件数字是巨大的。你的行列式是否也接近于 0？很有可能）。底线可能是这种矩阵的负幂对舍入和其他类型的错误非常敏感。除非在单独处理病态案例的任何一个程序中都有一些潜在的启发式。