ruby-changes:20305
From: marcandre <ko1@a...>
Date: Fri, 1 Jul 2011 15:30:29 +0900 (JST)
Subject: [ruby-changes:20305] marcandRe: r32353 (trunk): * lib/matrix: Add LUP decomposition
marcandre 2011-07-01 15:23:12 +0900 (Fri, 01 Jul 2011) New Revision: 32353 http://svn.ruby-lang.org/cgi-bin/viewvc.cgi?view=rev&revision=32353 Log: * lib/matrix: Add LUP decomposition Added files: trunk/lib/matrix/lup_decomposition.rb Modified files: trunk/ChangeLog trunk/NEWS trunk/lib/matrix.rb Index: ChangeLog =================================================================== --- ChangeLog (revision 32352) +++ ChangeLog (revision 32353) @@ -1,3 +1,7 @@ +Fri Jul 1 15:23:00 2011 Marc-Andre Lafortune <ruby-core@m...> + + * lib/matrix: Add LUP decomposition + Fri Jul 1 15:21:14 2011 Marc-Andre Lafortune <ruby-core@m...> * lib/matrix.rb: Allow non integer exponents for Matrix#** Index: lib/matrix/lup_decomposition.rb =================================================================== --- lib/matrix/lup_decomposition.rb (revision 0) +++ lib/matrix/lup_decomposition.rb (revision 32353) @@ -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 Index: lib/matrix.rb =================================================================== --- lib/matrix.rb (revision 32352) +++ lib/matrix.rb (revision 32353) @@ -98,6 +98,8 @@ # Matrix decompositions: # * <tt> #eigen </tt> # * <tt> #eigensystem </tt> +# * <tt> #lup </tt> +# * <tt> #lup_decomposition </tt> # # Complex arithmetic: # * <tt> conj </tt> @@ -122,6 +124,7 @@ include Enumerable include ExceptionForMatrix autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition" + autoload :LUPDecomposition, "matrix/lup_decomposition" # instance creations private_class_method :new @@ -1187,6 +1190,21 @@ 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 -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= #++ Index: NEWS =================================================================== --- NEWS (revision 32352) +++ NEWS (revision 32353) @@ -159,12 +159,15 @@ * 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? -- ML: ruby-changes@q... Info: http://www.atdot.net/~ko1/quickml/