非满秩矩阵也能求逆

今天遇到一个很奇怪的问题：一个方阵，逆矩阵存在，但不是满秩。

• 问题来源

    在实际应用的时候，发现返回值都是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的 源代码，发现引起这个问题的原因可能是精度问题，测试之后果不其然。

•  结论

    判断矩阵的逆矩阵是否存在时，一定要特别小心用满秩作为条件来判断，很可能会由于精度原因导致不可预估的结果。