Matlab mldivide转ArrayFire solve结果不一致问题求助

问题背景 #

我正在将Matlab代码迁移至C++的ArrayFire库，其中Matlab核心代码为：

tmp = (R\v0);

R是9×9复双精度矩阵，v0是9×1复双精度向量。ArrayFire中对应的实现为：

af::array tmp = af::solve(R, v0_ndx, AF_MAT_NONE);

（v0_ndx为ArrayFire中与v0完全一致的复向量），但两者输出结果差异显著，需排查原因并修正。

输入数据 #

R矩阵 #

(1.0000,0.0000)   (0.4330,0.5572)   (-0.4989,0.3674)  (0.6645,-0.1705)  (-0.5022,-0.3560) (0.0136,-0.5150)  (0.4447,0.2106)   (-0.5152,0.1762)  (0.7409,0.2369) 
(0.4330,-0.5572)  (1.0000,0.0000)   (0.2444,0.6848)   (0.3885,-0.7665)  (-0.5994,0.3623)  (-0.4474,-0.4617) (0.5153,-0.2913)  (-0.1954,0.6379)  (0.7954,-0.4984) 
(-0.4989,-0.3674) (0.2444,-0.6848)  (1.0000,0.0000)   (-0.4839,-0.1905) (0.4046,0.7260)   (-0.2262,0.3428)  (-0.1181,-0.5652) (0.6709,0.2600)   (-0.1537,-0.6409) 
(0.6645,0.1705)   (0.3885,0.7665)   (-0.4839,0.1905)  (1.0000,0.0000)   (-0.4047,-0.5696) (0.3418,-0.7446)  (0.2774,0.3234)   (-0.6889,-0.1392) (0.7275,0.4112) 
(-0.5022,0.3560)  (-0.5994,-0.3623) (0.4046,-0.7260)  (-0.4047,0.5696)  (1.0000,0.0000)   (0.3853,0.5248)   (-0.5786,0.0954)  (0.5647,-0.4580)  (-0.6453,0.0559) 
(0.0136,0.5150)   (-0.4474,0.4617)  (-0.2262,-0.3428) (0.3418,0.7446)   (0.3853,-0.5248)  (1.0000,0.0000)   (-0.2308,0.1934)  (-0.0831,-0.7353) (-0.1197,0.6637) 
(0.4447,-0.2106)  (0.5153,0.2913)   (-0.1181,0.5652)  (0.2774,-0.3234)  (-0.5786,-0.0954) (-0.2308,-0.1934) (1.0000,0.0000)   (-0.0868,0.3612)  (0.5367,-0.0571) 
(-0.5152,-0.1762) (-0.1954,-0.6379) (0.6709,-0.2600)  (-0.6889,0.1392)  (0.5647,0.4580)   (-0.0831,0.7353)  (-0.0868,-0.3612) (1.0000,0.0000)   (-0.4853,-0.3956) 
(0.7409,-0.2369)  (0.7954,0.4984)   (-0.1537,0.6409)  (0.7275,-0.4112)  (-0.6453,-0.0559) (-0.1197,-0.6637) (0.5367,0.0571)   (-0.4853,0.3956)  (1.0000,0.0000) 

v0/v0_ndx向量 #

(1.0000,0.0000) 
(0.8827,0.4699) 
(0.1784,-0.9840) 
(-0.9616,0.2745) 
(-0.5396,0.8419) 
(0.9908,-0.1353) 
(-0.5675,-0.8234) 
(0.8795,-0.4758) 
(0.9408,0.3389) 

运算结果对比 #

Matlab mldivide 结果 #

6.83220363990185 - 10.9619257689701i
40.3356738405241 + 25.0998371535142i
4.83385229841040 - 37.2045139160788i
-44.4504672789456 + 4.68752510309329i
-18.7242610414281 + 10.0126185742854i
13.2133070644273 - 1.29385372118540i
-2.14180173734413 + 12.1355411649223i
-4.16564154531280 - 6.04196278680202i
-2.43069655994337 + 6.59335067976014i

ArrayFire solve（AF_MAT_NONE）结果 #

(6.1867,-9.3346) 
(31.5626,-5.3650) 
(16.7719,-11.4021) 
(-31.0820,-10.4292) 
(-20.8629,35.2594) 
(3.3563,-6.7437) 
(-21.9321,-8.4480) 
(3.5929,5.8619) 
(-2.2422,9.7816) 

原因分析 #

1. 矩阵病态性：先计算R的条件数，若条件数极大（比如远大于1e10），说明矩阵是病态的，不同求解器的数值稳定性差异会导致结果出现显著偏差。Matlab的\运算符会自动检测矩阵状态并切换到更稳定的算法，而ArrayFire默认的LU分解对病态矩阵鲁棒性较差。
2. 求解算法差异：Matlab的\会根据矩阵的特性（是否对称、正定、病态等）自动选择LU、QR或Cholesky等分解算法；而ArrayFire的af::solve(..., AF_MAT_NONE)默认使用带部分选主元的LU分解，没有自动切换逻辑，当矩阵病态时容易出现数值误差。
3. 数据类型验证：需确认ArrayFire中R和v0_ndx是否为af::cf64（复双精度）类型，若误用单精度会导致精度损失。

修正方案 #

1. 改用QR分解求解：手动调用ArrayFire的QR分解，模拟Matlab对病态矩阵的处理逻辑：

   af::array Q, R_qr;
   af::qr(Q, R_qr, R); // 对R做QR分解
   af::array tmp = af::solve(R_qr, af::matmulTN(Q, v0_ndx), AF_MAT_UPPER); // 求解上三角矩阵系统
   

2. 矩阵正则化：若矩阵确实病态，可添加小的对角扰动增强稳定性：

   double alpha = 1e-6; // 可根据条件数调整
   af::array R_reg = R + alpha * af::eye(9, 9, af::cf64);
   af::array tmp = af::solve(R_reg, v0_ndx, AF_MAT_NONE);
   

3. 使用af::mldivide函数：如果你的ArrayFire版本支持（v3.8及以上），af::mldivide更接近Matlab的\行为，会自动选择最优求解算法：

   af::array tmp = af::mldivide(R, v0_ndx);
   

4. 验证数据类型：确保所有变量都使用复双精度：

   af::array R = af::array(9, 9, af::cf64, r_data); // r_data为复双精度数据数组
   af::array v0_ndx = af::array(9, 1, af::cf64, v_data);
   

内容的提问来源于stack exchange，提问作者Frank Allen