# Class: Matrix::LUPDecomposition

Inherits:
Object
• Object
show all
Includes:
ConversionHelper
Defined in:
lib/matrix/lup_decomposition.rb

## Overview

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.

## Instance Attribute Summary collapse

Returns the pivoting indices.

## Instance Method Summary collapse

• #det ⇒ Object (also: #determinant)

Returns the determinant of `A`, calculated efficiently from the factorization.

• constructor

A new instance of LUPDecomposition.

• Returns the permutation matrix `P`.

• Returns `true` if `U`, and hence `A`, is singular.

• Returns `m` so that `A*m = b`, or equivalently so that `L*U*m = P*b` `b` can be a Matrix or a Vector.

• #to_ary ⇒ Object (also: #to_a)

Returns `L`, `U`, `P` in an array.

• Returns the upper triangular factor `U`.

## Constructor Details

### #initialize(a) ⇒ LUPDecomposition

Returns a new instance of LUPDecomposition.

Raises:

• (TypeError)
 ``` 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217``` ```# File 'lib/matrix/lup_decomposition.rb', line 154 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_count = a.row_count @column_count = a.column_count @pivots = Array.new(@row_count) @row_count.times do |i| @pivots[i] = i end @pivot_sign = 1 lu_col_j = Array.new(@row_count) # Outer loop. @column_count.times do |j| # Make a copy of the j-th column to localize references. @row_count.times do |i| lu_col_j[i] = @lu[i][j] end # Apply previous transformations. @row_count.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_count-1) do |i| if (lu_col_j[i].abs > lu_col_j[p].abs) p = i end end if (p != j) @column_count.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_count && @lu[j][j] != 0) (j+1).upto(@row_count-1) do |i| @lu[i][j] = @lu[i][j].quo(@lu[j][j]) end end end end```

## Instance Attribute Details

Returns the pivoting indices

 ``` 63 64 65``` ```# File 'lib/matrix/lup_decomposition.rb', line 63 def pivots @pivots end```

## Instance Method Details

### #det ⇒ ObjectAlso known as: determinant

Returns the determinant of `A`, calculated efficiently from the factorization.

 ``` 79 80 81 82 83 84 85 86 87 88``` ```# File 'lib/matrix/lup_decomposition.rb', line 79 def det if (@row_count != @column_count) raise Matrix::ErrDimensionMismatch end d = @pivot_sign @column_count.times do |j| d *= @lu[j][j] end d end```

### #l ⇒ Object

 ``` 22 23 24 25 26 27 28 29 30 31 32``` ```# File 'lib/matrix/lup_decomposition.rb', line 22 def l Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j| if (i > j) @lu[i][j] elsif (i == j) 1 else 0 end end end```

### #p ⇒ Object

Returns the permutation matrix `P`

 ``` 48 49 50 51 52``` ```# File 'lib/matrix/lup_decomposition.rb', line 48 def p rows = Array.new(@row_count){Array.new(@row_count, 0)} @pivots.each_with_index{|p, i| rows[i][p] = 1} Matrix.send :new, rows, @row_count end```

### #singular? ⇒ Boolean

Returns `true` if `U`, and hence `A`, is singular.

Returns:

• (Boolean)
 ``` 67 68 69 70 71 72 73 74``` ```# File 'lib/matrix/lup_decomposition.rb', line 67 def singular? @column_count.times do |j| if (@lu[j][j] == 0) return true end end false end```

### #solve(b) ⇒ Object

Returns `m` so that `A*m = b`, or equivalently so that `L*U*m = P*b` `b` can be a Matrix or a Vector

 ``` 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152``` ```# File 'lib/matrix/lup_decomposition.rb', line 95 def solve b if (singular?) raise Matrix::ErrNotRegular, "Matrix is singular." end if b.is_a? Matrix if (b.row_count != @row_count) raise Matrix::ErrDimensionMismatch end # Copy right hand side with pivoting nx = b.column_count m = @pivots.map{|row| b.row(row).to_a} # Solve L*Y = P*b @column_count.times do |k| (k+1).upto(@column_count-1) do |i| nx.times do |j| m[i][j] -= m[k][j]*@lu[i][k] end end end # Solve U*m = Y (@column_count-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_count) raise Matrix::ErrDimensionMismatch end # Copy right hand side with pivoting m = b.values_at(*@pivots) # Solve L*Y = P*b @column_count.times do |k| (k+1).upto(@column_count-1) do |i| m[i] -= m[k]*@lu[i][k] end end # Solve U*m = Y (@column_count-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```

### #to_ary ⇒ ObjectAlso known as: to_a

Returns `L`, `U`, `P` in an array

 ``` 56 57 58``` ```# File 'lib/matrix/lup_decomposition.rb', line 56 def to_ary [l, u, p] end```

### #u ⇒ Object

Returns the upper triangular factor `U`

 ``` 36 37 38 39 40 41 42 43 44``` ```# File 'lib/matrix/lup_decomposition.rb', line 36 def u Matrix.build([@column_count, @row_count].min, @column_count) do |i, j| if (i <= j) @lu[i][j] else 0 end end end```