Skip to content
Browse files

* lib/matrix: Add LUP decomposition

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@32353 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
  • Loading branch information...
1 parent 59a3d59 commit c8e23882574a7cf14c5800369ca1e838eea7fe88 @marcandre marcandre committed Jul 1, 2011
Showing with 243 additions and 0 deletions.
  1. +4 −0 ChangeLog
  2. +3 −0 NEWS
  3. +18 −0 lib/matrix.rb
  4. +218 −0 lib/matrix/lup_decomposition.rb
View
4 ChangeLog
@@ -1,3 +1,7 @@
+Fri Jul 1 15:23:00 2011 Marc-Andre Lafortune <ruby-core@marc-andre.ca>
+
+ * lib/matrix: Add LUP decomposition
+
Fri Jul 1 15:21:14 2011 Marc-Andre Lafortune <ruby-core@marc-andre.ca>
* lib/matrix.rb: Allow non integer exponents for Matrix#**
View
3 NEWS
@@ -159,12 +159,15 @@ with all sufficient information, see the ChangeLog file.
* matrix
* new classes:
* Matrix::EigenvalueDecomposition
+ * Matrix::LUPDecomposition
* new methods:
* Matrix#diagonal?
* Matrix#eigen
* Matrix#eigensystem
* Matrix#hermitian?
* Matrix#lower_triangular?
+ * Matrix#lup
+ * Matrix#lup_decomposition
* Matrix#normal?
* Matrix#orthogonal?
* Matrix#permutation?
View
18 lib/matrix.rb
@@ -98,6 +98,8 @@ module ExceptionForMatrix # :nodoc:
# Matrix decompositions:
# * <tt> #eigen </tt>
# * <tt> #eigensystem </tt>
+# * <tt> #lup </tt>
+# * <tt> #lup_decomposition </tt>
#
# Complex arithmetic:
# * <tt> conj </tt>
@@ -122,6 +124,7 @@ class Matrix
include Enumerable
include ExceptionForMatrix
autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition"
+ autoload :LUPDecomposition, "matrix/lup_decomposition"
# instance creations
private_class_method :new
@@ -1187,6 +1190,21 @@ def eigensystem
end
alias eigen eigensystem
+ #
+ # Returns the LUP decomposition of the matrix; see +LUPDecomposition+.
+ # a = Matrix[[1, 2], [3, 4]]
+ # l, u, p = a.lup
+ # l.lower_triangular? # => true
+ # u.upper_triangular? # => true
+ # p.permutation? # => true
+ # l * u == a * p # => true
+ # a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)]
+ #
+ def lup
+ LUPDecomposition.new(self)
+ end
+ alias lup_decomposition lup
+
#--
# COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#++
View
218 lib/matrix/lup_decomposition.rb
@@ -0,0 +1,218 @@
+class Matrix
+ # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
+
+ #
+ # For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
+ # unit lower triangular matrix L, an n-by-n upper triangular matrix U,
+ # and a m-by-m permutation matrix P so that L*U = P*A.
+ # If m < n, then L is m-by-m and U is m-by-n.
+ #
+ # The LUP decomposition with pivoting always exists, even if the matrix is
+ # singular, so the constructor will never fail. The primary use of the
+ # LU decomposition is in the solution of square systems of simultaneous
+ # linear equations. This will fail if singular? returns true.
+ #
+
+ class LUPDecomposition
+ # Returns the lower triangular factor +L+
+
+ include Matrix::ConversionHelper
+
+ def l
+ Matrix.build(@row_size, @col_size) do |i, j|
+ if (i > j)
+ @lu[i][j]
+ elsif (i == j)
+ 1
+ else
+ 0
+ end
+ end
+ end
+
+ # Returns the upper triangular factor +U+
+
+ def u
+ Matrix.build(@col_size, @col_size) do |i, j|
+ if (i <= j)
+ @lu[i][j]
+ else
+ 0
+ end
+ end
+ end
+
+ # Returns the permutation matrix +P+
+
+ def p
+ rows = Array.new(@row_size){Array.new(@col_size, 0)}
+ @pivots.each_with_index{|p, i| rows[i][p] = 1}
+ Matrix.send :new, rows, @col_size
+ end
+
+ # Returns +L+, +U+, +P+ in an array
+
+ def to_ary
+ [l, u, p]
+ end
+ alias_method :to_a, :to_ary
+
+ # Returns the pivoting indices
+
+ attr_reader :pivots
+
+ # Returns +true+ if +U+, and hence +A+, is singular.
+
+ def singular? ()
+ @col_size.times do |j|
+ if (@lu[j][j] == 0)
+ return true
+ end
+ end
+ false
+ end
+
+ # Returns the determinant of +A+, calculated efficiently
+ # from the factorization.
+
+ def det
+ if (@row_size != @col_size)
+ Matrix.Raise Matrix::ErrDimensionMismatch unless square?
+ end
+ d = @pivot_sign
+ @col_size.times do |j|
+ d *= @lu[j][j]
+ end
+ d
+ end
+ alias_method :determinant, :det
+
+ # Returns +m+ so that <tt>A*m = b</tt>,
+ # or equivalently so that <tt>L*U*m = P*b</tt>
+ # +b+ can be a Matrix or a Vector
+
+ def solve b
+ if (singular?)
+ Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
+ end
+ if b.is_a? Matrix
+ if (b.row_size != @row_size)
+ Matrix.Raise Matrix::ErrDimensionMismatch
+ end
+
+ # Copy right hand side with pivoting
+ nx = b.column_size
+ m = @pivots.map{|row| b.row(row).to_a}
+
+ # Solve L*Y = P*b
+ @col_size.times do |k|
+ (k+1).upto(@col_size-1) do |i|
+ nx.times do |j|
+ m[i][j] -= m[k][j]*@lu[i][k]
+ end
+ end
+ end
+ # Solve U*m = Y
+ (@col_size-1).downto(0) do |k|
+ nx.times do |j|
+ m[k][j] = m[k][j].quo(@lu[k][k])
+ end
+ k.times do |i|
+ nx.times do |j|
+ m[i][j] -= m[k][j]*@lu[i][k]
+ end
+ end
+ end
+ Matrix.send :new, m, nx
+ else # same algorithm, specialized for simpler case of a vector
+ b = convert_to_array(b)
+ if (b.size != @row_size)
+ Matrix.Raise Matrix::ErrDimensionMismatch
+ end
+
+ # Copy right hand side with pivoting
+ m = b.values_at(*@pivots)
+
+ # Solve L*Y = P*b
+ @col_size.times do |k|
+ (k+1).upto(@col_size-1) do |i|
+ m[i] -= m[k]*@lu[i][k]
+ end
+ end
+ # Solve U*m = Y
+ (@col_size-1).downto(0) do |k|
+ m[k] = m[k].quo(@lu[k][k])
+ k.times do |i|
+ m[i] -= m[k]*@lu[i][k]
+ end
+ end
+ Vector.elements(m, false)
+ end
+ end
+
+ def initialize a
+ raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
+ # Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+ @lu = a.to_a
+ @row_size = a.row_size
+ @col_size = a.column_size
+ @pivots = Array.new(@row_size)
+ @row_size.times do |i|
+ @pivots[i] = i
+ end
+ @pivot_sign = 1
+ lu_col_j = Array.new(@row_size)
+
+ # Outer loop.
+
+ @col_size.times do |j|
+
+ # Make a copy of the j-th column to localize references.
+
+ @row_size.times do |i|
+ lu_col_j[i] = @lu[i][j]
+ end
+
+ # Apply previous transformations.
+
+ @row_size.times do |i|
+ lu_row_i = @lu[i]
+
+ # Most of the time is spent in the following dot product.
+
+ kmax = [i, j].min
+ s = 0
+ kmax.times do |k|
+ s += lu_row_i[k]*lu_col_j[k]
+ end
+
+ lu_row_i[j] = lu_col_j[i] -= s
+ end
+
+ # Find pivot and exchange if necessary.
+
+ p = j
+ (j+1).upto(@row_size-1) do |i|
+ if (lu_col_j[i].abs > lu_col_j[p].abs)
+ p = i
+ end
+ end
+ if (p != j)
+ @col_size.times do |k|
+ t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
+ end
+ k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
+ @pivot_sign = -@pivot_sign
+ end
+
+ # Compute multipliers.
+
+ if (j < @row_size && @lu[j][j] != 0)
+ (j+1).upto(@row_size-1) do |i|
+ @lu[i][j] = @lu[i][j].quo(@lu[j][j])
+ end
+ end
+ end
+ end
+ end
+end

0 comments on commit c8e2388

Please sign in to comment.
Something went wrong with that request. Please try again.