[sdpb] 30/233: Added cholesky update functions
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:14 UTC 2017
This is an automated email from the git hooks/post-receive script.
thansen pushed a commit to branch master
in repository sdpb.
commit 375a288268db611b99952f38e55e16e0851fb2a0
Author: David Simmons-Duffin <dsd at athena.sns.ias.edu>
Date: Thu Jul 17 15:59:15 2014 -0400
Added cholesky update functions
---
Makefile | 2 +-
main.cpp | 157 ++++++++++++++++++++++++++++++++++++------------------
mpack/Rrot.cpp | 99 ++++++++++++++++++++++++++++++++++
mpack/Rrotg.cpp | 103 +++++++++++++++++++++++++++++++++++
mpack/mblas_gmp.h | 2 +
5 files changed, 310 insertions(+), 53 deletions(-)
diff --git a/Makefile b/Makefile
index 7412c03..cf36984 100755
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-OBJECTS = main.o tinyxml2.o mpack/Rdot.o \
+OBJECTS = main.o tinyxml2.o mpack/Rrotg.o mpack/Rrot.o mpack/Rdot.o \
mpack/Rtrsm.o mpack/Rsyrk.o mpack/Raxpy.o \
mpack/Rgemm.o mpack/Rtrmm.o mpack/Rtrsv.o \
mpack/iMlaenv.o mpack/Rlamch.o mpack/Rlascl.o \
diff --git a/main.cpp b/main.cpp
index 7e4675e..7d622c1 100644
--- a/main.cpp
+++ b/main.cpp
@@ -19,8 +19,6 @@ using std::ofstream;
using tinyxml2::XMLDocument;
using tinyxml2::XMLElement;
-Real Infinity = 1e200;
-
template <class T>
ostream& operator<<(ostream& os, const vector<T>& v) {
os << "{";
@@ -301,72 +299,89 @@ Real frobeniusProductOfSums(const Matrix &X, const Matrix &dX,
return result;
}
-
-// Not currently supporting this. Should probably switch to mpfr...
-//
-// void matrixMultiplyFirstSym(Matrix &A, Matrix &B, Matrix &result) {
-// assert(A.cols == A.rows);
-// assert(A.cols == B.rows);
-// assert(B.rows == result.rows);
-// assert(B.cols == result.cols);
-
-// Rsymm("Left", "Upper", B.rows, B.cols, 1,
-// &A.elements[0], A.rows,
-// &B.elements[0], B.rows,
-// 0,
-// &result.elements[0], result.rows);
-// }
-// result = choleskyDecomposition(a) (lower triangular)
+// L (lower triangular) such that A = L L^T
// Inputs:
-// - a : dim x dim symmetric matrix
-// - result : dim x dim lower-triangular matrix
+// - A : dim x dim symmetric matrix
+// - L : dim x dim lower-triangular matrix
//
-void choleskyDecomposition(Matrix &a, Matrix &result) {
- int dim = a.rows;
- assert(a.cols == dim);
- assert(result.rows == dim);
- assert(result.cols == dim);
+void choleskyDecomposition(Matrix &A, Matrix &L) {
+ int dim = A.rows;
+ assert(A.cols == dim);
+ assert(L.rows == dim);
+ assert(L.cols == dim);
+ // Set lower-triangular part of L to cholesky decomposition
+ L.copyFrom(A);
mpackint info;
- Real *resultArray = &result.elements[0];
-
- Rcopy(dim*dim, &a.elements[0], 1, resultArray, 1);
-
- // The lower-triangular part of result is now our cholesky matrix
- Rpotrf("Lower", dim, resultArray, dim, &info);
+ Rpotrf("Lower", dim, &L.elements[0], dim, &info);
assert(info == 0);
- // Set the upper-triangular part of the result to zero
+ // Set the upper-triangular part of the L to zero
for (int j = 0; j < dim; j++)
for (int i = 0; i < j; i++)
- result.elements[i + j*dim] = 0;
+ L.elements[i + j*dim] = 0;
}
-// result = a^-1
+// L' (lower triangular) such that L' L'^T = L L^T + v v^T. i.e., if L
+// is a cholesky decomposition of A, then L' is a cholesky
+// decomposition of A + v v^T
+// Inputs:
+// - L : dim x dim lower-triangular matrix
+// - v : pointer to the head of a length-dim vector
+// both L and v are modified in place
+//
+void choleskyUpdate(Matrix &L, Real *v) {
+ int dim = L.rows;
+ Real c, s, x, y;
+ for (int r = 0; r < dim; r++) {
+ x = L.get(r,r);
+ y = *(v+r);
+ Rrotg(&x, &y, &c, &s);
+ Rrot(dim - r, &L.elements[r*(dim+1)], 1, v+r, 1, c, s);
+ }
+}
+
+// L' (lower triangular) such that L' L'^T = L L^T + V V^T. i.e., if L
+// is a cholesky decomposition of A, then L' is a cholesky
+// decomposition of A + V V^T. This is more efficient than directly
+// computing the cholesky decomposition of A + V V^T if V has a small
+// number of columns.
+// Inputs:
+// - L : dim x dim lower-triangular matrix
+// - V : dim x n matrix
+// both L and V are modified in place
+//
+void choleskyUpdate(Matrix &L, Matrix &V) {
+ assert(L.rows == V.rows);
+ for (int c = 0; c < V.cols; c++)
+ choleskyUpdate(L, &V.elements[c*V.rows]);
+}
+
+// result = A^-1
// Inputs:
-// - a : dim x dim lower-triangular matrix
+// - A : dim x dim lower-triangular matrix
// - result : dim x dim lower-triangular matrix
//
-void inverseLowerTriangular(Matrix &a, Matrix &result) {
- int dim = a.rows;
- assert(a.cols == dim);
+void inverseLowerTriangular(Matrix &A, Matrix &result) {
+ int dim = A.rows;
+ assert(A.cols == dim);
assert(result.rows == dim);
assert(result.cols == dim);
result.setIdentity();
Rtrsm("Left", "Lower", "NoTranspose", "NonUnitDiagonal",
- dim, dim, 1, &a.elements[0], dim, &result.elements[0], dim);
+ dim, dim, 1, &A.elements[0], dim, &result.elements[0], dim);
}
-// result = choleskyDecomposition(a)^-1
+// result = choleskyDecomposition(A)^-1
// Inputs:
-// - a : dim x dim symmetric matrix
+// - A : dim x dim symmetric matrix
// - work : dim x dim matrix
// - result : dim x dim lower-triangular matrix
//
-void inverseCholesky(Matrix &a, Matrix &work, Matrix &result) {
- choleskyDecomposition(a, work);
+void inverseCholesky(Matrix &A, Matrix &work, Matrix &result) {
+ choleskyDecomposition(A, work);
inverseLowerTriangular(work, result);
}
@@ -1051,13 +1066,13 @@ public:
Real lambdaStar;
int maxIterations;
SolverParameters():
- betaStar(0.1),
- betaBar(0.3),
- epsilonStar(1e-30),
- epsilonDash(1e-30),
- gammaStar(0.9),
- alphaMax(100),
- lambdaStar(1e4),
+ betaStar("0.1"),
+ betaBar("0.3"),
+ epsilonStar("1e-30"),
+ epsilonDash("1e-30"),
+ gammaStar("0.9"),
+ alphaMax("100"),
+ lambdaStar("1e4"),
maxIterations(100) {}
};
@@ -1468,7 +1483,7 @@ Real minEigenvalueViaQR(BlockDiagonalMatrix &A, vector<Vector> &eigenvalues, vec
assert(A.blocks.size() == eigenvalues.size());
assert(A.blocks.size() == workspace.size());
- Real lambdaMin = Infinity;
+ Real lambdaMin = A.diagonalPart[0];
for (unsigned int i = 0; i < A.diagonalPart.size(); i++)
lambdaMin = min(lambdaMin, A.diagonalPart[i]);
@@ -1894,6 +1909,43 @@ void testSDPSolver(const char *file) {
datfile.close();
}
+void testCholeskyUpdate() {
+ Matrix A(4,4);
+ Matrix B(A);
+ Matrix C(A);
+ Matrix L(A);
+ Matrix LT(L);
+ Matrix V(4, 2);
+ Matrix VT(V.cols, V.rows);
+ V.set(0,0,1);
+ V.set(1,0,2);
+ V.set(2,0,3);
+ V.set(3,0,4);
+ V.set(0,1,5);
+ V.set(1,1,4);
+ V.set(2,1,3);
+ V.set(3,1,2);
+ for (int r = 0; r < V.rows; r++)
+ for (int c = 0; c < V.cols; c++)
+ VT.set(c, r, V.get(r,c));
+ Matrix U(V);
+
+ A.addDiagonal(4);
+ cout << "A = " << A << endl;
+ cout << "V = " << V << endl;
+ choleskyDecomposition(A, L);
+ choleskyUpdate(L, U);
+ LT = L;
+ LT.transpose();
+
+ matrixMultiply(V, VT, B);
+ B += A;
+ matrixMultiply(L, LT, C);
+ C -= B;
+
+ cout << "L L^T - (A + V V^T) = " << C << endl;
+}
+
int main(int argc, char** argv) {
mpf_set_default_prec(200);
@@ -1902,7 +1954,8 @@ int main(int argc, char** argv) {
//testBlockCongruence();
//testBlockDiagonalCholesky();
- testSDPSolver(argv[1]);
+ //testSDPSolver(argv[1]);
+ testCholeskyUpdate();
//testMinEigenvalue();
diff --git a/mpack/Rrot.cpp b/mpack/Rrot.cpp
new file mode 100644
index 0000000..c51e28f
--- /dev/null
+++ b/mpack/Rrot.cpp
@@ -0,0 +1,99 @@
+/*
+ * Copyright (c) 2008-2010
+ * Nakata, Maho
+ * All rights reserved.
+ *
+ * $Id: Rrot.cpp,v 1.5 2010/08/07 05:50:10 nakatamaho Exp $
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ */
+/*
+Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
+ *
+ * $Id: Rrot.cpp,v 1.5 2010/08/07 05:50:10 nakatamaho Exp $
+
+$COPYRIGHT$
+
+Additional copyrights may follow
+
+$HEADER$
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+Based on http://www.netlib.org/blas/drot.f
+applies a plane rotation.
+*/
+
+#include <mblas_gmp.h>
+
+void Rrot(mpackint n, mpf_class * dx, mpackint incx, mpf_class * dy, mpackint incy, mpf_class c, mpf_class s)
+{
+ mpackint i, ix, iy;
+ mpf_class temp;
+
+ if (n <= 0)
+ return;
+ ix = 0;
+ iy = 0;
+
+ if (incx < 0)
+ ix = (-n + 1) * incx;
+ if (incy < 0)
+ iy = (-n + 1) * incy;
+ for (i = 0; i < n; i++) {
+ temp = c * dx[ix] + s * dy[iy];
+ dy[iy] = c * dy[iy] - s * dx[ix];
+ dx[ix] = temp;
+ ix = ix + incx;
+ iy = iy + incy;
+ }
+ return;
+}
diff --git a/mpack/Rrotg.cpp b/mpack/Rrotg.cpp
new file mode 100644
index 0000000..b310497
--- /dev/null
+++ b/mpack/Rrotg.cpp
@@ -0,0 +1,103 @@
+/*
+ * Copyright (c) 2008-2010
+ * Nakata, Maho
+ * All rights reserved.
+ *
+ * $Id: Rrotg.cpp,v 1.5 2010/08/07 05:50:10 nakatamaho Exp $
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ */
+/*
+Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
+ *
+ * $Id: Rrotg.cpp,v 1.5 2010/08/07 05:50:10 nakatamaho Exp $
+
+$COPYRIGHT$
+
+Additional copyrights may follow
+
+$HEADER$
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+Based on http://www.netlib.org/blas/drotg.f
+*/
+
+#include <mblas_gmp.h>
+
+void Rrotg(mpf_class * da, mpf_class * db, mpf_class * c, mpf_class * s)
+{
+ mpf_class roe, scale;
+ mpf_class r, z;
+
+ roe = *db;
+ if (abs(*da) > abs(*db))
+ roe = *da;
+ scale = abs(*da) + abs(*db);
+ if (scale == 0.0) {
+ (*c) = 1.0;
+ (*s) = 0.0;
+ r = 0.0;
+ z = 0.0;
+ } else {
+ r = scale * sqrt(((*da) / scale) * ((*da) / scale) + ((*db) / scale) * ((*db) / scale));
+ r = Msign(1.0, roe) * r;
+ (*c) = (*da) / r;
+ (*s) = (*db) / r;
+ z = 1.0;
+ if (abs(*da) > abs(*db))
+ z = *s;
+ if (abs(*db) >= abs(*da) && *c != 0.0)
+ z = 1.0 / (*c);
+ }
+ *da = r;
+ *db = z;
+}
diff --git a/mpack/mblas_gmp.h b/mpack/mblas_gmp.h
index 87863de..0c0320b 100644
--- a/mpack/mblas_gmp.h
+++ b/mpack/mblas_gmp.h
@@ -46,6 +46,8 @@
_MPACK_EXTERN_ int mpack_errno;
/* LEVEL 1 MBLAS */
+void Rrotg(mpf_class * da, mpf_class * db, mpf_class * c, mpf_class * s);
+void Rrot(mpackint n, mpf_class * dx, mpackint incx, mpf_class * dy, mpackint incy, mpf_class c, mpf_class s);
mpf_class Rdot(mpackint n, mpf_class * dx, mpackint incx, mpf_class * dy,
mpackint incy);
void Rcopy(mpackint n, mpf_class * dx, mpackint incx, mpf_class * dy,
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/sdpb.git
More information about the debian-science-commits
mailing list