1. 程式人生 > >非滿秩矩陣也能求逆

非滿秩矩陣也能求逆

今天遇到一個很奇怪的問題:一個方陣,逆矩陣存在,但不是滿秩。

  • 問題來源

    在實際應用的時候,發現返回值都是0,於是跟蹤到這裡,發現了這個問題:JtJ不是滿秩,因此JtJN保持初始化的零值。

mat JtJN = zeros(N, N);
mat Result = zeros(N, 1);
if (N == rank(JtJ))
{
    JtJN = inv(JtJ);
}
for (int i = 0; i < N; i++)
{
    Result[i] = 1.96 * sqrt(re * abs(JtJN(i, i)) / M) / (p0[i] + EPS);
}
return Result;
  • 問題跟蹤

    我把JtJ這個資料儲存成txt,然後在Matlab裡load進去。

1.42564385120145	46664.0071835617	-361.723628266118	1.43453776778475	2.27890158958207	1491.23699249997	2.66016443927362
46664.0071835617	2561658858.23584	-30028903.1298037	47411.4041209533	107097.080643255	129526675.413329	143227.334019439
-361.723628266118	-30028903.1298037	464526.466675988	-373.437398450597	-1274.63678807808	-2315792.08994344	-1969.46141027209
1.43453776778475	47411.4041209533	-373.437398450597	1.44372953400922	2.31104869171119	1551.42174553270	2.71062044898402
2.27890158958207	107097.080643255	-1274.63678807808	2.31104869171119	5.20067695242223	7813.59509020671	7.98532576281474
1491.23699249997	129526675.413329	-2315792.08994344	1551.42174553270	7813.59509020671	25339933.2510429	20655.0780322129
2.66016443927362	143227.334019439	-1969.46141027209	2.71062044898402	7.98532576281474	20655.0780322129	50.0000010000000

    請看以下結果:

    於是,我跟蹤進rank的 原始碼,發現引起這個問題的原因可能是精度問題,測試之後果不其然。

  •  結論

    判斷矩陣的逆矩陣是否存在時,一定要特別小心用滿秩作為條件來判斷,很可能會由於精度原因導致不可預估的結果。