numpy.linalg.solve でのエラー対処法について

Question

import numpy as np
A = np.array(
[[147, 420, 21]
 [420, 1202, 60]
 [21, 60, 3]])
b = np.array([-9443, -27060, -1349])
d = np.linalg.det(A)
x = np.linalg.solve(A, b)

Aの行列式は本来0ですが，
d = 1.30562227696e-13
と出力されてしまいます.おかげでsolveはnp.linalg.linalg.LinAlgErrorをraiseしてくれず，とんちんかんな答えを出します.どうしたらよいでしょうか？

Answer 1

値が1.30562227696e-13ととんでもなく0に近いので、どこかで計算誤差が生じています。

numpyのリファレンスを見てみましょう。
numpy.linalg.detには、次のような記述があります。

See also:
slogdet Another way to representing the determinant, more suitable for large matrices where underflow/overflow may occur.

このアドバイスに従って、numpy.linalg.slogdetを使ってみます。

>>> np.linalg.slogdet(A)
(1.0, -30.984227930851471)

リファレンスを見ると、この結果の捉え方がわかります。

In all cases, the determinant is equal to sign * np.exp(logdet).

つまり、行列式はnumpy.e**-31と計算されたわけです。
slogdetの力を借りても、正しく計算することは出来ませんでした。

---

でも、悲観する必要はありません。
解が不定であることを確かめるためには、階数を計算すればよいです。

>>> np.linalg.matrix_rank(A)
2

solveを呼び出す前に階数を計算して分岐させましょう。

---

あるいは、次のような強引な方法も案外柔軟かもしれません。

try:
    if np.linalg.matrix_rank(A) < A.shape[0]:
        raise np.linalg.linalg.LinAlgError
    else:
        np.linalg.solve(A, b)
except:
    print('解は不定です')

Comments

[hayatoU]

rankのことを忘れていました。助かりました。