e$B$^$D$b$He(B e$B$f$-$R$m$G$9e(B
rubyforgee$B%P%0%l%]!<%H$K0J2<$N$h$&$J%l%]!<%H$,Mh$F$^$9$,!";de(B
e$B$K$OH=CG$G$-$^$;$s!#$I$&$G$7$g$&!#@PDM$5$s!)e(B
e$B%3%a%s%H$OJT=8$7$F$"$j$^$9!#%Q%C%A$Oe(B(3e$B$D$"$ke(B)e$B:G=*HG$N$b$Ne(B
e$B$G$9!#e(B
There is a problem with the Matrix.inverse method.
Given the following matrix
m = Matrix[[0.417787968829298, 0.542037605627772,
0.729142268138943],
[0.0, 6.12303176911189e-17, 1.0],
[-0.542037605627772, 0.417787968829298, -2.55812900589452e-17]]
the result of
m.inverse * m
is
[[1.0, 0.0, 0.650423523448542],
[5.55111512312578e-17, 1.0, 0.84385869292008],
[0.0, 9.24446373305873e-33, 1.0]]
but it should be
[[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0]]
or close to this!
The patch below should “fix” the inverse method. It adds partial
pivoting to the Gauss-Jordan algorithm, making it stable. (I used Ruby
1.8.5 (2006-08-22).)
I don’t think there is an error in the algorithm, and the matrix is not
ill-conditioned (a condition number of about 2.5 isn’t so bad).
Looks to me like BigDecimal uses Gaussian elimination with partial
pivoting too, or at least that’s what the comments say. It includes some
additional scaling, but this may be specific to BigDecimal (precision
works differently for BigDecimals and floats).
Peter
— lib/matrix.rb.old 2007-03-19 01:46:44.000000000 +0100
+++ lib/matrix.rb 2007-03-19 03:26:38.000000000 +0100
@@ -590,15 +590,21 @@
a = src.to_a
for k in 0..size
-
if (akk = a[k][k]) == 0 -
i = k -
begin -
Matrix.Raise ErrNotRegular if (i += 1) > size -
end while a[i][k] == 0
-
i = k -
akk = a[k][k].abs -
for j in (k+1)..size -
v = a[j][k].abs -
if v > akk -
i = j -
akk = v -
end -
end -
Matrix.Raise ErrNotRegular if akk == 0 -
if i != k a[i], a[k] = a[k], a[i] @rows[i], @rows[k] = @rows[k], @rows[i]
-
akk = a[k][k] end
-
akk = a[k][k] for i in 0 .. size next if i == k