LU decomposition: Difference between revisions
Line 571: | Line 571: | ||
=={{header|J}}== |
=={{header|J}}== |
||
Taken with slight modification from [http://www.jsoftware.com/jwiki/Essays/LU%20Decomposition] |
Taken with slight modification from [http://www.jsoftware.com/jwiki/Essays/LU%20Decomposition]. |
||
'''Solution:''' |
'''Solution:''' |
||
<lang j> mp=: +/ .* |
<lang j> mp=: +/ .* |
Revision as of 01:28, 21 June 2011
You are encouraged to solve this task according to the task description, using any language you may know.
Every square matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle A} can be decomposed into a product of a lower triangular matrix and a upper triangular matrix , as described in LU decomposition.
It is a modified form of Gaussian elimination. While the Cholesky decomposition only works for symmetric, positive definite matrices, the more general LU decomposition works for any square matrix.
There are several algorithms for calculating L and U. To derive Crout's algorithm for a 3x3 example, we have to solve the following system:
We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of are set to 1
so we get a solvable system of 9 unknowns and 9 equations.
Solving for the other and , we get the following equations:
and for :
We see that there is a calculation pattern, which can be expressed as the following formulas, first for
and then for
We see in the second formula that to get the below the diagonal, we have to divide by the diagonal element (pivot) , so we get problems when is either 0 or very small, which leads to numerical instability.
The solution to this problem is pivoting , which means rearranging the rows of , prior to the decomposition, in a way that the largest element of each column gets onto the diagonal of . Rearranging the columns means to multiply by a permutation matrix :
Example:
The decomposition algorithm is then applied on the rearranged matrix so that
Task description
The task is to implement a routine which will take a square nxn matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle A} and return a lower triangular matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle L} , a upper triangular matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle U} and a permutation matrix Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle P} , so that the above equation is fullfilled. You should then test it on the following two examples and include your output.
Example 1:
A 1 3 5 2 4 7 1 1 0 L 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 U 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 P 0 1 0 1 0 0 0 0 1
Example 2:
A 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 L 1.00000 0.00000 0.00000 0.00000 0.27273 1.00000 0.00000 0.00000 0.09091 0.28750 1.00000 0.00000 0.18182 0.23125 0.00360 1.00000 U 11.00000 9.00000 24.00000 2.00000 0.00000 14.54545 11.45455 0.45455 0.00000 0.00000 -3.47500 5.68750 0.00000 0.00000 0.00000 0.51079 P 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
Ada
decomposition.ads: <lang Ada>with Ada.Numerics.Generic_Real_Arrays; generic
with package Matrix is new Ada.Numerics.Generic_Real_Arrays (<>);
package Decomposition is
-- decompose a square matrix A by PA = LU procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix);
end Decomposition;</lang>
decomposition.adb: <lang Ada>package body Decomposition is
procedure Swap_Rows (M : in out Matrix.Real_Matrix; From, To : Natural) is Temporary : Matrix.Real; begin if From = To then return; end if; for I in M'Range (2) loop Temporary := M (M'First (1) + From, I); M (M'First (1) + From, I) := M (M'First (1) + To, I); M (M'First (1) + To, I) := Temporary; end loop; end Swap_Rows;
function Pivoting_Matrix (M : Matrix.Real_Matrix) return Matrix.Real_Matrix is use type Matrix.Real; Order : constant Positive := M'Length (1); Result : Matrix.Real_Matrix := Matrix.Unit_Matrix (Order); Max : Matrix.Real; Row : Natural; begin for J in 0 .. Order - 1 loop Max := M (M'First (1) + J, M'First (2) + J); Row := J; for I in J .. Order - 1 loop if M (M'First (1) + I, M'First (2) + J) > Max then Max := M (M'First (1) + I, M'First (2) + J); Row := I; end if; end loop; if J /= Row then -- swap rows J and Row Swap_Rows (Result, J, Row); end if; end loop; return Result; end Pivoting_Matrix;
procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix) is use type Matrix.Real_Matrix, Matrix.Real; Order : constant Positive := A'Length (1); A2 : Matrix.Real_Matrix (A'Range (1), A'Range (2)); S : Matrix.Real; begin L := (others => (others => 0.0)); U := (others => (others => 0.0)); P := Pivoting_Matrix (A); A2 := P * A; for J in 0 .. Order - 1 loop L (L'First (1) + J, L'First (2) + J) := 1.0; for I in 0 .. J loop S := 0.0; for K in 0 .. I - 1 loop S := S + U (U'First (1) + K, U'First (2) + J) * L (L'First (1) + I, L'First (2) + K); end loop; U (U'First (1) + I, U'First (2) + J) := A2 (A2'First (1) + I, A2'First (2) + J) - S; end loop; for I in J + 1 .. Order - 1 loop S := 0.0; for K in 0 .. J loop S := S + U (U'First (1) + K, U'First (2) + J) * L (L'First (1) + I, L'First (2) + K); end loop; L (L'First (1) + I, L'First (2) + J) := (A2 (A2'First (1) + I, A2'First (2) + J) - S) / U (U'First (1) + J, U'First (2) + J); end loop; end loop; end Decompose;
end Decomposition;</lang>
Example usage: <lang Ada>with Ada.Numerics.Real_Arrays; with Ada.Text_IO; with Decomposition; procedure Decompose_Example is
package Real_Decomposition is new Decomposition (Matrix => Ada.Numerics.Real_Arrays);
package Real_IO is new Ada.Text_IO.Float_IO (Float); procedure Print (M : Ada.Numerics.Real_Arrays.Real_Matrix) is begin for Row in M'Range (1) loop for Col in M'Range (2) loop Real_IO.Put (M (Row, Col), 3, 2, 0); end loop; Ada.Text_IO.New_Line; end loop; end Print;
Example_1 : constant Ada.Numerics.Real_Arrays.Real_Matrix := ((1.0, 3.0, 5.0), (2.0, 4.0, 7.0), (1.0, 1.0, 0.0)); P_1, L_1, U_1 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_1'Range (1), Example_1'Range (2)); Example_2 : constant Ada.Numerics.Real_Arrays.Real_Matrix := ((11.0, 9.0, 24.0, 2.0), (1.0, 5.0, 2.0, 6.0), (3.0, 17.0, 18.0, 1.0), (2.0, 5.0, 7.0, 1.0)); P_2, L_2, U_2 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_2'Range (1), Example_2'Range (2));
begin
Real_Decomposition.Decompose (A => Example_1, P => P_1, L => L_1, U => U_1); Real_Decomposition.Decompose (A => Example_2, P => P_2, L => L_2, U => U_2); Ada.Text_IO.Put_Line ("Example 1:"); Ada.Text_IO.Put_Line ("A:"); Print (Example_1); Ada.Text_IO.Put_Line ("L:"); Print (L_1); Ada.Text_IO.Put_Line ("U:"); Print (U_1); Ada.Text_IO.Put_Line ("P:"); Print (P_1); Ada.Text_IO.New_Line; Ada.Text_IO.Put_Line ("Example 2:"); Ada.Text_IO.Put_Line ("A:"); Print (Example_2); Ada.Text_IO.Put_Line ("L:"); Print (L_2); Ada.Text_IO.Put_Line ("U:"); Print (U_2); Ada.Text_IO.Put_Line ("P:"); Print (P_2);
end Decompose_Example;</lang>
Output:
Example 1: A: 1.00 3.00 5.00 2.00 4.00 7.00 1.00 1.00 0.00 L: 1.00 0.00 0.00 0.50 1.00 0.00 0.50 -1.00 1.00 U: 2.00 4.00 7.00 0.00 1.00 1.50 0.00 0.00 -2.00 P: 0.00 1.00 0.00 1.00 0.00 0.00 0.00 0.00 1.00 Example 2: A: 11.00 9.00 24.00 2.00 1.00 5.00 2.00 6.00 3.00 17.00 18.00 1.00 2.00 5.00 7.00 1.00 L: 1.00 0.00 0.00 0.00 0.27 1.00 0.00 0.00 0.09 0.29 1.00 0.00 0.18 0.23 0.00 1.00 U: 11.00 9.00 24.00 2.00 0.00 14.55 11.45 0.45 0.00 0.00 -3.47 5.69 0.00 0.00 0.00 0.51 P: 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00
Common Lisp
Uses the routine (mmul A B) from Matrix multiplication.
<lang lisp>;; Creates a nxn identity matrix. (defun eye (n)
(let ((I (make-array `(,n ,n) :initial-element 0))) (loop for j from 0 to (- n 1) do (setf (aref I j j) 1)) I))
- Swap two rows l and k of a mxn matrix A, which is a 2D array.
(defun swap-rows (A l k)
(let* ((n (cadr (array-dimensions A))) (row (make-array n :initial-element 0))) (loop for j from 0 to (- n 1) do (setf (aref row j) (aref A l j)) (setf (aref A l j) (aref A k j)) (setf (aref A k j) (aref row j)))))
- Creates the pivoting matrix for A.
(defun pivotize (A)
(let* ((n (car (array-dimensions A))) (P (eye n))) (loop for j from 0 to (- n 1) do (let ((max (aref A j j)) (row j)) (loop for i from j to (- n 1) do (if (> (aref A i j) max) (setq max (aref A i j) row i))) (if (not (= j row)) (swap-rows P j row))))
;; Return P. P))
- Decomposes a square matrix A by PA=LU and returns L, U and P.
(defun lu (A)
(let* ((n (car (array-dimensions A))) (L (make-array `(,n ,n) :initial-element 0)) (U (make-array `(,n ,n) :initial-element 0)) (P (pivotize A)) (A (mmul P A)))
(loop for j from 0 to (- n 1) do (setf (aref L j j) 1) (loop for i from 0 to j do (setf (aref U i j) (- (aref A i j) (loop for k from 0 to (- i 1) sum (* (aref U k j) (aref L i k)))))) (loop for i from j to (- n 1) do (setf (aref L i j) (/ (- (aref A i j) (loop for k from 0 to (- j 1) sum (* (aref U k j) (aref L i k)))) (aref U j j)))))
;; Return L, U and P. (values L U P)))</lang>
Example 1:
<lang lisp>(setf g (make-array '(3 3) :initial-contents '((1 3 5) (2 4 7)(1 1 0))))
- 2A((1 3 5) (2 4 7) (1 1 0))
(lu g)
- 2A((1 0 0) (1/2 1 0) (1/2 -1 1))
- 2A((2 4 7) (0 1 3/2) (0 0 -2))
- 2A((0 1 0) (1 0 0) (0 0 1))</lang>
Example 2:
<lang lisp>(setf h (make-array '(4 4) :initial-contents '((11 9 24 2)(1 5 2 6)(3 17 18 1)(2 5 7 1))))
- 2A((11 9 24 2) (1 5 2 6) (3 17 18 1) (2 5 7 1))
(lup h)
- 2A((1 0 0 0) (3/11 1 0 0) (1/11 23/80 1 0) (2/11 37/160 1/278 1))
- 2A((11 9 24 2) (0 160/11 126/11 5/11) (0 0 -139/40 91/16) (0 0 0 71/139))
- 2A((1 0 0 0) (0 0 1 0) (0 1 0 0) (0 0 0 1))</lang>
D
Using some functions from Matrix multiplication task.
<lang d>import std.stdio, std.algorithm, std.typecons;
/// Creates a nxn identity matrix. pure T[][] identityMatrix(T)(const int n)
in { assert(n >= 0); } body { auto m = new typeof(return)(n, n); foreach (i, row; m) { row[] = 0; row[i] = 1; } return m; }
/// Creates the pivoting matrix for m. pure T[][] pivotize(T)(const T[][] m)
in { // is square assert(isRectangular(m) && m[0].length == m.length); } body { int n = m.length; auto P = identityMatrix!T(n);
foreach (j; 0 .. n) { T max = m[j][j]; int row = j; foreach (i; j .. n) if (m[i][j] > max) { max = m[i][j]; row = i; } if (j != row) swap(P[j], P[row]); }
return P; }
/// Decomposes a square matrix A by PA=LU and returns L, U and P. /*pure*/ Tuple!(T[][], T[][], T[][]) lu(T)(const T[][] A)
in { // is square assert(isRectangular(A) && A[0].length == A.length); } body { int n = A.length; auto L = new T[][](n, n); auto U = new T[][](n, n); foreach (i; 0 .. n) { L[i][i .. $] = 0; U[i][0 .. i] = 0; } auto P = pivotize!T(A); auto A2 = matrixMul!T(P, A);
foreach (j; 0 .. n) { L[j][j] = 1; foreach (i; 0 .. j+1) { T s1 = 0; foreach (k; 0 .. i) s1 += U[k][j] * L[i][k]; U[i][j] = A2[i][j] - s1; } foreach (i; j .. n) { T s2 = 0; foreach (k; 0 .. j) s2 += U[k][j] * L[i][k]; L[i][j] = (A2[i][j] - s2) / U[j][j]; } }
return tuple(L, U, P); }
void main() {
enum double[][] a = [[1, 3, 5], [2, 4, 7], [1, 1, 0]]; foreach (part; lu(a).tupleof) writeln(prettyPrint(part), "\n"); writeln();
enum double[][] b = [[11, 9, 24, 2], [1, 5, 2, 6], [3, 17, 18, 1], [2, 5, 7, 1]]; foreach (part; lu(b).tupleof) writeln(prettyPrint(part), "\n");
}</lang> Output:
[[1, 0, 0], [0.5, 1, 0], [0.5, -1, 1]] [[2, 4, 7], [0, 1, 1.5], [0, 0, -2]] [[0, 1, 0], [1, 0, 0], [0, 0, 1]] [[1, 0, 0, 0], [0.272727, 1, 0, 0], [0.090909, 0.2875, 1, 0], [0.181818, 0.23125, 0.00359712, 1]] [[11, 9, 24, 2], [0, 14.5455, 11.4545, 0.454545], [0, 0, -3.475, 5.6875], [0, 0, 0, 0.510791]] [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]
J
Taken with slight modification from [1].
Solution: <lang j> mp=: +/ .*
LU=: 3 : 0
'm n'=. $ A=. y if. 1=m do.
p ; (=1) ; p{"1 A [ p=. C. (n-1);~.0,(0~:,A)i.1
else.
m2=. >.m%2 'p1 L1 U1'=. LU m2{.A D=. (/:p1) {"1 m2}.A F=. m2 {."1 D E=. m2 {."1 U1 FE1=. F mp %. E G=. m2}."1 D - FE1 mp U1 'p2 L2 U2'=. LU G p3=. (i.m2),m2+p2 H=. (/:p3) {"1 U1 (p1{p3) ; (L1,FE1,.L2) ; H,(-n){."1 U2
end. )
permtomat=. 1 {.~"0 -@>:@:/: LUdecompose=: (permtomat&.>@{. , }.)@:LU</lang>
Example use: <lang j> A=.3 3$1 3 5 2 4 7 1 1 0
LUdecompose A
┌─────┬─────┬───────┐ │1 0 0│1 0 0│1 3 5│ │0 1 0│2 1 0│0 _2 _3│ │0 0 1│1 1 1│0 0 _2│ └─────┴─────┴───────┘
mp/> LUdecompose A
1 3 5 2 4 7 1 1 0</lang>
Python
<lang python>from pprint import pprint
def matrixMul(A, B):
TB = zip(*B) return [[sum(ea*eb for ea,eb in zip(a,b)) for b in TB] for a in A]
def pivotize(m):
"""Creates the pivoting matrix for m.""" n = len(m) ID = [[float(i == j) for i in xrange(n)] for j in xrange(n)] for j in xrange(n): row = max(xrange(j, n), key=lambda i: m[i][j]) if j != row: ID[j], ID[row] = ID[row], ID[j] return ID
def lu(A):
"""Decomposes a nxn matrix A by PA=LU and returns L, U and P.""" n = len(A) L = [[0.0] * n for i in xrange(n)] U = [[0.0] * n for i in xrange(n)] P = pivotize(A) A2 = matrixMul(P, A) for j in xrange(n): L[j][j] = 1.0 for i in xrange(j+1): s1 = sum(U[k][j] * L[i][k] for k in xrange(i)) U[i][j] = A2[i][j] - s1 for i in xrange(j, n): s2 = sum(U[k][j] * L[i][k] for k in xrange(j)) L[i][j] = (A2[i][j] - s2) / U[j][j] return (L, U, P)
a = [[1, 3, 5], [2, 4, 7], [1, 1, 0]] for part in lu(a):
pprint(part, width=19) print
print b = [[11,9,24,2],[1,5,2,6],[3,17,18,1],[2,5,7,1]] for part in lu(b):
pprint(part) print</lang>
Output:
[[1.0, 0.0, 0.0], [0.5, 1.0, 0.0], [0.5, -1.0, 1.0]] [[2.0, 4.0, 7.0], [0.0, 1.0, 1.5], [0.0, 0.0, -2.0]] [[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]] [[1.0, 0.0, 0.0, 0.0], [0.27272727272727271, 1.0, 0.0, 0.0], [0.090909090909090912, 0.28749999999999998, 1.0, 0.0], [0.18181818181818182, 0.23124999999999996, 0.0035971223021580693, 1.0]] [[11.0, 9.0, 24.0, 2.0], [0.0, 14.545454545454547, 11.454545454545455, 0.45454545454545459], [0.0, 0.0, -3.4749999999999996, 5.6875], [0.0, 0.0, 0.0, 0.51079136690647597]] [[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
Tcl
<lang tcl>package require Tcl 8.5 namespace eval matrix {
namespace path {::tcl::mathfunc ::tcl::mathop}
# Construct an identity matrix of the given size proc identity {order} {
set m [lrepeat $order [lrepeat $order 0]] for {set i 0} {$i < $order} {incr i} { lset m $i $i 1 } return $m
}
# Produce the pivot matrix for a given matrix proc pivotize {matrix} {
set n [llength $matrix] set p [identity $n] for {set j 0} {$j < $n} {incr j} { set max [lindex $matrix $j $j] set row $j for {set i $j} {$i < $n} {incr i} { if {[lindex $matrix $i $j] > $max} { set max [lindex $matrix $i $j] set row $i } } if {$j != $row} { # Row swap inlined; too trivial to have separate procedure set tmp [lindex $p $j] lset p $j [lindex $p $row] lset p $row $tmp } } return $p
}
# Decompose a square matrix A by PA=LU and return L, U and P proc luDecompose {A} {
set n [llength $A] set L [lrepeat $n [lrepeat $n 0]] set U $L set P [pivotize $A] set A [multiply $P $A]
for {set j 0} {$j < $n} {incr j} { lset L $j $j 1 for {set i 0} {$i <= $j} {incr i} { lset U $i $j [- [lindex $A $i $j] [SumMul $L $U $i $j $i]] } for {set i $j} {$i < $n} {incr i} { set sum [SumMul $L $U $i $j $j] lset L $i $j [/ [- [lindex $A $i $j] $sum] [lindex $U $j $j]] } }
return [list $L $U $P]
}
# Helper that makes inner loop nicer; multiplies column and row, # possibly partially... proc SumMul {A B i j kmax} {
set s 0.0 for {set k 0} {$k < $kmax} {incr k} { set s [+ $s [* [lindex $A $i $k] [lindex $B $k $j]]] } return $s
}
}</lang> Support code: <lang tcl># Code adapted from Matrix_multiplication and Matrix_transposition tasks namespace eval matrix {
# Get the size of a matrix; assumes that all rows are the same length, which # is a basic well-formed-ness condition... proc size {m} {
set rows [llength $m] set cols [llength [lindex $m 0]] return [list $rows $cols]
}
# Matrix multiplication implementation proc multiply {a b} {
lassign [size $a] a_rows a_cols lassign [size $b] b_rows b_cols if {$a_cols != $b_rows} { error "incompatible sizes: a($a_rows, $a_cols), b($b_rows, $b_cols)" } set temp [lrepeat $a_rows [lrepeat $b_cols 0]] for {set i 0} {$i < $a_rows} {incr i} { for {set j 0} {$j < $b_cols} {incr j} { lset temp $i $j [SumMul $a $b $i $j $a_cols] } } return $temp
}
# Pretty printer for matrices proc print {matrix {fmt "%g"}} {
set max [Widest $matrix $fmt] lassign [size $matrix] rows cols foreach row $matrix { foreach val $row width $max { puts -nonewline [format "%*s " $width [format $fmt $val]] } puts "" }
} proc Widest {m fmt} {
lassign [size $m] rows cols set max [lrepeat $cols 0] foreach row $m { for {set j 0} {$j < $cols} {incr j} { set s [format $fmt [lindex $row $j]] lset max $j [max [lindex $max $j] [string length $s]] } } return $max
}
}</lang> Demonstrating: <lang tcl># This does the decomposition and prints it out nicely proc demo {A} {
lassign $A L U P foreach v {A L U P} {
upvar 0 $v matrix puts "${v}:" matrix::print $matrix %.5g if {$v ne "P"} {puts "---------------------------------"}
}
} demo {{1 3 5} {2 4 7} {1 1 0}} puts "=================================" demo {{11 9 24 2} {1 5 2 6} {3 17 18 1} {2 5 7 1}}</lang> Output:
A: 1 3 5 2 4 7 1 1 0 --------------------------------- L: 1 0 0 0.5 1 0 0.5 -1 1 --------------------------------- U: 2 4 7 0 1 1.5 0 0 -2 --------------------------------- P: 0 1 0 1 0 0 0 0 1 ================================= A: 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 --------------------------------- L: 1 0 0 0 0.27273 1 0 0 0.090909 0.2875 1 0 0.18182 0.23125 0.0035971 1 --------------------------------- U: 11 9 24 2 0 14.545 11.455 0.45455 0 0 -3.475 5.6875 0 0 0 0.51079 --------------------------------- P: 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1