From 17329b05128f48d1c30f6278540cc2a2e8c56c14 Mon Sep 17 00:00:00 2001 From: Matthew Peddie Date: Mon, 30 May 2016 19:41:51 +1000 Subject: [PATCH 1/2] EXAMPLES: SPD matrix decompositions - LDL^T decomposition from Golub and Van Loan - UDU^T decomposition from Gibbs - Outer-product Cholesky decomposition from Golub and Van Loan - Rank-1 Cholesky update/downdate from Golub and Van Loan --- examples/cholesky.plv | 103 ++++++++++++++++++ examples/cholesky_test.plv | 209 +++++++++++++++++++++++++++++++++++++ 2 files changed, 312 insertions(+) create mode 100644 examples/cholesky.plv create mode 100644 examples/cholesky_test.plv diff --git a/examples/cholesky.plv b/examples/cholesky.plv new file mode 100644 index 0000000..73f9fdd --- /dev/null +++ b/examples/cholesky.plv @@ -0,0 +1,103 @@ +import prelude; + +-- "Matrix Computations" 4th ed., Golub and Van Loan. Algorithm +-- 4.1.2, the LDL^T decomposition for a symmetric, positive-definite +-- matrix. +-- +-- On success, the argument A now contains the lower-triangular factor +-- L and the diagonal vector D. +ldlt {m} + (inout A :: double[m, m]) + :: s8 := ( + + v :: double[m]; + + for j in 0:m -> ( + -- Compute v[1:j]. + for i in 0:j -> ( + v[i] <- A[j, i] * A[i, i]; + ); + v[j] <- A[j, j] - A[j, 0:j] * v[0:j]; + + -- Store d[j]. + A[j, j] <- v[j]; + + -- Compute L[j+1:n, j]. + A[j+1:m, j] <- (A[j+1:m, j] - A[j+1:m, 0:j]*v[0:j]) / v[j]; + ); + + return 0; +); + +-- "Advanced Kalman Filtering, Least-Squares and Modeling.", Bruce +-- P. Gibbs. Algorithm 10.2-2, the UDU^T decomposition for a +-- symmetric, positive-definite matrix. +-- +-- On exit, the argument A now contains the upper-triangular factor U +-- and the diagonal vector D. +udut {m} + (inout A :: double[m, m]) + :: Void := ( + for j in m-1:0:-1 -> ( -- column loop + alpha := if (A[j, j] > 0.0) then (1.0 / A[j, j]) else 0.0; + for k in 0:j -> ( -- row loop + beta := A[k, j]; + A[k, j] <- alpha*beta; -- scale elements in j column + A[0..k, k] <- A[0..k, k] - beta*A[0..k, j]; -- modify column k + ); + ); +); + + +-- "Matrix Computations" 4th ed., Golub and Van Loan. Algorithm +-- 4.2.2, the Cholesky decomposition for a symmetric, +-- positive-definite matrix. +-- +-- On exit, the argument A now contains the lower-triangular +-- Cholesky factor G. +cholesky {m} + (inout A :: double[m, m]) + :: Void := ( + for k in 0:m -> ( + sq := sqrt(A[k, k]); + A[k, k] <- sq; + A[k+1:m, k] <- A[k+1:m, k] / sq; + for j in k+1:m -> ( + A[j:m, j] <- A[j:m, j] - A[j:m, k]*A[j, k]; + ); + ); +); + +-- "Matrix computations" 4th ed., Golub and Van Loan. Based on +-- Algorithm 4.2.2, the Cholesky decomposition by outer product +-- updates. +-- +-- On exit, the lower-triangular Cholesky factor G passed as an +-- argument has been updated to the lower-triangular Cholesky factor +-- G' of A + vv^T, where A = GG^T. +-- +-- If a negative 'sign' argument is given, then the rank-1 Cholesky +-- downdate A - vv^T is performed instead. It is the caller's +-- responsibility to ensure that GG^T remains positive definite. +cholupdate {m} + (inout G :: double[m, m]) + (v' :: double[m]) + (sign :: s8) + :: Void := ( + + v := v'; + + for k in 0:m -> ( + sq := sqrt(G[k, k]^2 + (if sign < 0 + then -v[k]^2 + else v[k]^2)); + alpha := sq / G[k, k]; + beta := v[k] / G[k, k]; + + G[k, k] <- sq; + G[k+1:m, k] <- (G[k+1:m, k] + (if sign < 0 + then -beta*v[k+1:m] + else beta*v[k+1:m])) / alpha; + v[k+1:m] <- alpha*v[k+1:m] - beta*G[k+1:m, k]; + ); +); diff --git a/examples/cholesky_test.plv b/examples/cholesky_test.plv new file mode 100644 index 0000000..b9e2ec6 --- /dev/null +++ b/examples/cholesky_test.plv @@ -0,0 +1,209 @@ +import cholesky; + +static diagonal {n} + (d :: double[n]) + (out m :: double[n, n]) + :: Void := ( + for i in 0:n -> + for j in 0:n -> + if i == j then (m[i, j] <- d[i]) else (m[i, j] <- 0.0); +); + +static tril {m, n} + (a :: double[m, n]) + (out b :: double[m, n]) + :: Void := ( + for i in 0:m -> + for j in 0:n -> + if i < j then b[i, j] <- 0.0 else b[i, j] <- a[i, j]; +); + +static triu {m, n} + (a :: double[m, n]) + (out b :: double[m, n]) + :: Void := ( + for i in 0:m -> + for j in 0:n -> + if i > j then b[i, j] <- 0.0 else b[i, j] <- a[i, j]; +); + +static udut_test () :: s32 := ( + u0 := mat(1.0, 0.57733, 0.42103, 0.98750; + 0.0, 1.0, 0.99883, 0.99510; + 0.0, 0.0, 1.0, 0.53597; + 0.0, 0.0, 0.00000, 1.0); + d0 := mat(2.0, 0, 0, 0; + 0, 3.0, 0, 0; + 0, 0, 4.0, 0; + 0, 0, 0, 5.0); + + A0 := u0*d0*u0^T; + + A1 := A0; + + udut(inout A1); + + d :: double[4]; + d <- vec j in 4 -> A1[j, j]; + u :: double[4, 4]; + triu A1 (out u); + for i in 0:4 -> u[i, i] <- 1.0; + + dd :: double[4, 4]; + diagonal d (out dd); + + printf "\nresults:\nA:\n"; + print_mat A0; + printf "u0:\n"; + print_mat u0; + printf "u:\n"; + print_mat u; + printf "d:\n"; + print_mat dd; + + derror := d0 - dd; + uerror := u0 - u; + tol := 1e-9; + for col in 2 -> ( + assert(norm(derror[:, col]) < tol); + assert(norm(uerror[:, col]) < tol); + ); + + printf "Solution to within %3e!\n" tol; +); + +static ldlt_test () :: s32 := ( + l0 := mat(1.0, 0.0, 0.0, 0.0; + 0.37526, 1.0, 0.0, 0.0; + 0.25545, 0.35379, 1.0, 0.0; + 0.24456, 0.79429, 0.44717, 1.0); + d0 := mat(5.0, 0, 0, 0; + 0, 4.0, 0, 0; + 0, 0, 3.0, 0; + 0, 0, 0, 2.0); + + A0 := l0*d0*l0^T; + + A1 := A0; + + ret := ldlt(inout A1); + + d :: double[4]; + d <- vec j in 4 -> A1[j, j]; + l :: double[4, 4]; + tril A1 (out l); + for i in 0:4 -> l[i, i] <- 1.0; + + dd :: double[4, 4]; + diagonal d (out dd); + + printf "\nresults:\nA:\n"; + print_mat A0; + printf "l0:\n"; + print_mat l0; + printf "l:\n"; + print_mat l; + printf "d:\n"; + print_mat dd; + + derror := d0 - dd; + lerror := l0 - l; + tol := 1e-9; + for col in 4 -> ( + assert(norm(derror[:, col]) < tol); + assert(norm(lerror[:, col]) < tol); + ); + + printf "Solution to within %3e!\n" tol; + + ret; +); + +static cholesky_test () :: s32 := ( + g0 := mat(1.58518, 0.00000, 0.00000, 0.00000; + 1.18082, 0.22361, 0.00000, 0.00000; + -1.22629, -1.09190, 1.10581, 0.00000; + -0.19586, -0.61511, 1.34394, 0.82566); + + A0 := g0*g0^T; + + A1 := A0; + + cholesky(inout A1); + + g :: double[4, 4]; + tril A1 (out g); + + printf "\nresults:\nA:\n"; + print_mat A0; + printf "g0:\n"; + print_mat g0; + printf "g:\n"; + print_mat g; + + gerror := g0 - g; + tol := 1e-9; + for col in 4 -> ( + assert(norm(gerror[:, col]) < tol); + ); + + printf "Solution to within %3e!\n" tol; +); + +static cholupdate_test () :: s32 := ( + g0 := mat(1.58518, 0.00000, 0.00000, 0.00000; + 1.18082, 0.22361, 0.00000, 0.00000; + -1.22629, -1.09190, 1.10581, 0.00000; + -0.19586, -0.61511, 1.34394, 0.82566); + + v0 := vec(0.33595, 0.24187, 0.36134, 0.42910); + + A0 := g0*g0^T; + + gup' := g0; + gdown' := g0; + + cholupdate (inout gup') v0 1; + cholupdate (inout gdown') v0 (-1); + + Aup := A0 + v0*v0^T; + cholesky(inout Aup); + gup :: double[4, 4]; + tril Aup (out gup); + + Adown := A0 - v0*v0^T; + cholesky(inout Adown); + gdown :: double[4, 4]; + tril Adown (out gdown); + + printf "\nresults:\ng0:\n"; + print_mat g0; + printf "gup:\n"; + print_mat gup; + printf "gup':\n"; + print_mat gup'; + printf "gdown:\n"; + print_mat gdown; + printf "gdown':\n"; + print_mat gdown'; + + uperror := gup - gup'; + downerror := gdown - gdown'; + tol := 1e-9; + for col in 4 -> ( + assert(norm(uperror[:, col]) < tol); + assert(norm(downerror[:, col]) < tol); + ); + + printf "Solution to within %3e!\n" tol; + +); + +main () :: s32 := ( + ret := ldlt_test (); + udut_test (); + cholesky_test (); + cholupdate_test (); + + ret; +); \ No newline at end of file From c4da23df51e67cbe4c0a8dca14d587831e549fe0 Mon Sep 17 00:00:00 2001 From: Matthew Peddie Date: Mon, 30 May 2016 21:05:32 +1000 Subject: [PATCH 2/2] Add Cholesky decomposition tests. --- test/Main.hs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Main.hs b/test/Main.hs index f89d22b..a63a0a7 100644 --- a/test/Main.hs +++ b/test/Main.hs @@ -44,6 +44,7 @@ doCase (label, args, fn, outcome) = testCase label $ do cases :: [TestCase] cases = [ ("qr", ["-I", "examples"], "examples/qr_test.plv", (True, Just ExitSuccess)) + , ("cholesky", ["-I", "examples"], "examples/cholesky_test.plv", (True, Just ExitSuccess)) , ("cyclic", ["-I", "examples/module-tests"], "examples/module-tests/cycleA.plv", (False, Nothing)) ] @@ -52,4 +53,3 @@ tests :: TestTree tests = testGroup "Tests" [ testGroup "Units" $ map doCase cases ] -