1berry/ruby
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

* lib/matrix.rb: Improve algorithm for Matrix#determinant and Matrix#rank

Browse Source 
{determinant,det,rank}_e are now deprecated. [ruby-core:28273]
  Also fixes a bug in Determinant#rank (e.g. [[0,1][0,1][0,1]])

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@27554 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
This commit is contained in:
marcandre 2010-04-29 18:19:12 +00:00
parent a3a4542fb4
commit 0a3c78face
3 changed files with 118 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
ChangeLog
View File

@ -1,3 +1,12 @@
Fri Apr 30 03:17:20 2010 Marc-Andre Lafortune <ruby-core@marc-andre.ca>
* lib/matrix.rb: Improve algorithm for Matrix#determinant and
Matrix#rank
{determinant,det,rank}_e are now deprecated. [ruby-core:28273]
Also fixes a bug in Determinant#rank (e.g. [[0,1][0,1][0,1]])
Matrix#singular?, Matrix#regular? now raise on rectangular matrices
and use determinant instead of rank.
Fri Apr 30 00:52:56 2010 NAKAMURA Usaku <usa@ruby-lang.org> Fri Apr 30 00:52:56 2010 NAKAMURA Usaku <usa@ruby-lang.org>
* win32/Makefile.sub (config.h): define some constants to select * win32/Makefile.sub (config.h): define some constants to select

232
lib/matrix.rb
View File

@ -740,170 +740,146 @@ class Matrix
# #
# Returns the determinant of the matrix. # Returns the determinant of the matrix.
# This method's algorithm is Gaussian elimination method #
# and using Numeric#quo(). Beware that using Float values, with their # Beware that using Float values can yield erroneous results
# usual lack of precision, can affect the value returned by this method. Use # because of their lack of precision.
# Rational values or Matrix#det_e instead if this is important to you. # Consider using exact types like Rational or BigDecimal instead.
# #
# Matrix[[7,6], [3,9]].determinant # Matrix[[7,6], [3,9]].determinant
# => 45.0 # => 45
# #
def determinant def determinant
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
m = @rows
case row_size
# Up to 4x4, give result using Laplacian expansion by minors.
# This will typically be faster, as well as giving good results
# in case of Floats
when 0
+1
when 1
+ m[0][0]
when 2
+ m[0][0] * m[1][1] - m[0][1] * m[1][0]
when 3
m0 = m[0]; m1 = m[1]; m2 = m[2]
+ m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \
- m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \
+ m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
when 4
m0 = m[0]; m1 = m[1]; m2 = m[2]; m3 = m[3]
+ m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \
- m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \
+ m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \
- m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \
+ m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \
- m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \
+ m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \
- m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \
+ m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \
- m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \
+ m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \
- m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
else
# For bigger matrices, use an efficient and general algorithm.
# Currently, we use the Gauss-Bareiss algorithm
determinant_bareiss
end
end
alias_method :det, :determinant
#
# Private. Use Matrix#determinant
#
# Returns the determinant of the matrix, using
# Bareiss' multistep integer-preserving gaussian elimination.
# It has the same computational cost order O(n^3) as standard Gaussian elimination.
# Intermediate results are fraction free and of lower complexity.
# A matrix of Integers will have thus intermediate results that are also Integers,
# with smaller bignums (if any), while a matrix of Float will usually have more
# precise intermediate results.
#
def determinant_bareiss
size = row_size size = row_size
last = size - 1
a = to_a a = to_a
no_pivot = Proc.new{ return 0 }
det = 1 sign = +1
pivot = 1
size.times do |k| size.times do |k|
if (akk = a[k][k]) == 0 previous_pivot = pivot
i = (k+1 ... size).find {|ii| if (pivot = a[k][k]) == 0
a[ii][k] != 0 switch = (k+1 ... size).find(no_pivot) {|row|
a[row][k] != 0
} }
return 0 if i.nil? a[switch], a[k] = a[k], a[switch]
a[i], a[k] = a[k], a[i] pivot = a[k][k]
akk = a[k][k] sign = -sign
det *= -1
end end
(k+1).upto(last) do |i|
(k + 1 ... size).each do |ii| ai = a[i]
q = a[ii][k].quo(akk) (k+1).upto(last) do |j|
(k + 1 ... size).each do |j| ai[j] = (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
a[ii][j] -= a[k][j] * q
end end
end end
det *= akk
end end
det sign * pivot
end end
alias det determinant private :determinant_bareiss
# #
# Returns the determinant of the matrix. # deprecated; use Matrix#determinant
# This method's algorithm is Gaussian elimination method.
# This method uses Euclidean algorithm. If all elements are integer,
# really exact value. But, if an element is a float, can't return
# exact value.
#
# Matrix[[7,6], [3,9]].determinant
# => 63
# #
def determinant_e def determinant_e
Matrix.Raise ErrDimensionMismatch unless square? warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant"
rank
size = row_size
a = to_a
det = 1
size.times do |k|
if a[k][k].zero?
i = (k+1 ... size).find {|ii|
a[ii][k] != 0
}
return 0 if i.nil?
a[i], a[k] = a[k], a[i]
det *= -1
end
(k + 1 ... size).each do |ii|
q = a[ii][k].quo(a[k][k])
(k ... size).each do |j|
a[ii][j] -= a[k][j] * q
end
unless a[ii][k].zero?
a[ii], a[k] = a[k], a[ii]
det *= -1
redo
end
end
det *= a[k][k]
end
det
end end
alias det_e determinant_e alias det_e determinant_e
# #
# Returns the rank of the matrix. Beware that using Float values, # Returns the rank of the matrix.
# probably return faild value. Use Rational values or Matrix#rank_e # Beware that using Float values can yield erroneous results
# for getting exact result. # because of their lack of precision.
# Consider using exact types like Rational or BigDecimal instead.
# #
# Matrix[[7,6], [3,9]].rank # Matrix[[7,6], [3,9]].rank
# => 2 # => 2
# #
def rank def rank
if column_size > row_size # We currently use Bareiss' multistep integer-preserving gaussian elimination
a = transpose.to_a # (see comments on determinant)
a_column_size = row_size a = to_a
a_row_size = column_size last_column = column_size - 1
else last_row = row_size - 1
a = to_a
a_column_size = column_size
a_row_size = row_size
end
rank = 0 rank = 0
a_column_size.times do |k| pivot_row = 0
if (akk = a[k][k]) == 0 previous_pivot = 1
i = (k+1 ... a_row_size).find {|ii| 0.upto(last_column) do |k|
a[ii][k] != 0 switch_row = (pivot_row .. last_row).find {|row|
} a[row][k] != 0
if i }
a[i], a[k] = a[k], a[i] if switch_row
akk = a[k][k] a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
else pivot = a[pivot_row][k]
i = (k+1 ... a_column_size).find {|ii| (pivot_row+1).upto(last_row) do |i|
a[k][ii] != 0 ai = a[i]
} (k+1).upto(last_column) do |j|
next if i.nil? ai[j] = (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
(k ... a_column_size).each do |j| end
a[j][k], a[j][i] = a[j][i], a[j][k] end
end pivot_row += 1
akk = a[k][k] previous_pivot = pivot
end
end end
(k + 1 ... a_row_size).each do |ii|
q = a[ii][k].quo(akk)
(k + 1... a_column_size).each do |j|
a[ii][j] -= a[k][j] * q
end
end
rank += 1
end end
return rank pivot_row
end end
# #
# Returns the rank of the matrix. This method uses Euclidean # deprecated; use Matrix#rank
# algorithm. If all elements are integer, really exact value. But, if
# an element is a float, can't return exact value.
#
# Matrix[[7,6], [3,9]].rank
# => 2
# #
def rank_e def rank_e
a = to_a warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank"
a_column_size = column_size rank
a_row_size = row_size
pi = 0
a_column_size.times do |j|
if i = (pi ... a_row_size).find{|i0| !a[i0][j].zero?}
if i != pi
a[pi], a[i] = a[i], a[pi]
end
(pi + 1 ... a_row_size).each do |k|
q = a[k][j].quo(a[pi][j])
(pi ... a_column_size).each do |j0|
a[k][j0] -= q * a[pi][j0]
end
if k > pi && !a[k][j].zero?
a[k], a[pi] = a[pi], a[k]
redo
end
end
pi += 1
end
end
pi
end end

23
test/matrix/test_matrix.rb
View File

@ -250,14 +250,12 @@ class TestMatrix < Test::Unit::TestCase
assert(Matrix[[1, 0], [0, 1]].regular?) assert(Matrix[[1, 0], [0, 1]].regular?)
assert(Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].regular?) assert(Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].regular?)
assert(!Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].regular?) assert(!Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].regular?)
assert(!Matrix[[1, 0, 0], [0, 1, 0]].regular?)
end end
def test_singular? def test_singular?
assert(!Matrix[[1, 0], [0, 1]].singular?) assert(!Matrix[[1, 0], [0, 1]].singular?)
assert(!Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].singular?) assert(!Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].singular?)
assert(Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].singular?) assert(Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].singular?)
assert(Matrix[[1, 0, 0], [0, 1, 0]].singular?)
end end
def test_square? def test_square?
@ -325,31 +323,20 @@ class TestMatrix < Test::Unit::TestCase
end end
def test_det def test_det
assert_in_delta(45.0, Matrix[[7,6],[3,9]].det, 0.0001) assert_equal(45, Matrix[[7,6],[3,9]].det)
assert_in_delta(0.0, Matrix[[0,0],[0,0]].det, 0.0001) assert_equal(0, Matrix[[0,0],[0,0]].det)
assert_in_delta(-7.0, Matrix[[0,0,1],[0,7,6],[1,3,9]].det, 0.0001) assert_equal(-7, Matrix[[0,0,1],[0,7,6],[1,3,9]].det)
end assert_equal(42, Matrix[[7,0,1,0,12],[8,1,1,9,1],[4,0,0,-7,17],[-1,0,0,-4,8],[10,1,1,8,6]].det)
def test_det_e
assert_equal(45, Matrix[[7,6],[3,9]].det_e)
assert_equal(0, Matrix[[0,0],[0,0]].det_e)
assert_equal(-7, Matrix[[0,0,1],[0,7,6],[1,3,9]].det_e)
end end
def test_rank2 def test_rank2
assert_equal(2, Matrix[[7,6],[3,9]].rank) assert_equal(2, Matrix[[7,6],[3,9]].rank)
assert_equal(0, Matrix[[0,0],[0,0]].rank) assert_equal(0, Matrix[[0,0],[0,0]].rank)
assert_equal(3, Matrix[[0,0,1],[0,7,6],[1,3,9]].rank) assert_equal(3, Matrix[[0,0,1],[0,7,6],[1,3,9]].rank)
assert_equal(1, Matrix[[0,1],[0,1],[0,1]].rank)
assert_equal(2, @m1.rank) assert_equal(2, @m1.rank)
end end
def test_rank_e
assert_equal(2, Matrix[[7,6],[3,9]].rank_e)
assert_equal(0, Matrix[[0,0],[0,0]].rank_e)
assert_equal(3, Matrix[[0,0,1],[0,7,6],[1,3,9]].rank_e)
assert_equal(2, @m1.rank_e)
end
def test_trace def test_trace
assert_equal(1+5+9, Matrix[[1,2,3],[4,5,6],[7,8,9]].trace) assert_equal(1+5+9, Matrix[[1,2,3],[4,5,6],[7,8,9]].trace)
end end