[sdpb] 96/233: Now including free variables in the differential equation; way more numerically stable -- yay :)
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:24 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 9580423c8c436d640d90b37c19454e4e24fabc9e
Author: David Simmons-Duffin <dsd at athena.sns.ias.edu>
Date: Tue Oct 21 18:41:34 2014 -0400
Now including free variables in the differential equation; way more numerically stable -- yay :)
---
src/Matrix.cpp | 144 +-------------------------------------------
src/Matrix.h | 36 -----------
src/SDPSolver.cpp | 168 +++++++++++++++-------------------------------------
src/SDPSolver.h | 18 ++----
src/SDPSolverIO.cpp | 16 ++---
src/main.cpp | 6 +-
src/serialize.h | 3 +-
src/tests.cpp | 40 -------------
8 files changed, 68 insertions(+), 363 deletions(-)
diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index 72463da..0353bdd 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -47,24 +47,6 @@ void matrixMultiply(Matrix &A, Matrix &B, Matrix &C) {
matrixScaleMultiplyAdd(1, A, B, 0, C);
}
-// B = A^T A
-void matrixSquare(Matrix &A, Matrix &B) {
- assert(A.cols == B.rows);
- assert(B.cols == B.rows);
-
- #pragma omp parallel for schedule(dynamic)
- for (int c = 0; c < B.cols; c++) {
- for (int r = 0; r <= c; r++) {
- Real tmp = 0;
- for (int p = 0; p < A.rows; p++)
- tmp += A.elt(p,r) * A.elt(p,c);
- B.elt(r,c) = tmp;
- if (r != c)
- B.elt(c,r) = tmp;
- }
- }
-}
-
// Set block starting at (bRow, bCol) of B to A^T A
void matrixSquareIntoBlock(Matrix &A, Matrix &B, int bRow, int bCol) {
assert(bRow + A.cols <= B.rows);
@@ -83,22 +65,6 @@ void matrixSquareIntoBlock(Matrix &A, Matrix &B, int bRow, int bCol) {
}
}
-// A := L A L^T
-void lowerTriangularCongruence(Matrix &A, Matrix &L) {
- int dim = A.rows;
- assert(A.cols == dim);
- assert(L.rows == dim);
- assert(L.cols == dim);
-
- Rtrmm("Right", "Lower", "Transpose", "NonUnitDiagonal", dim, dim, 1,
- &L.elements[0], dim,
- &A.elements[0], dim);
-
- Rtrmm("Left", "Lower", "NoTranspose", "NonUnitDiagonal", dim, dim, 1,
- &L.elements[0], dim,
- &A.elements[0], dim);
-}
-
// A := L^{-1} A L^{-T}
void lowerTriangularInverseCongruence(Matrix &A, Matrix &L) {
int dim = A.rows;
@@ -217,22 +183,12 @@ void LUDecomposition(Matrix &A, vector<Integer> &pivots) {
assert(info == 0);
}
-void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Real *B, int bcols, int ldb) {
- Integer info;
- Rgetrs("NoTranspose", LU.rows, bcols, &LU.elements[0], LU.rows, &pivots[0], B, ldb, &info);
- assert(info == 0);
-}
-
-void solveWithLUDecompositionTranspose(Matrix &LU, vector<Integer> &pivots, Real *B, int bcols, int ldb) {
+void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b) {
Integer info;
- Rgetrs("Transpose", LU.rows, bcols, &LU.elements[0], LU.rows, &pivots[0], B, ldb, &info);
+ Rgetrs("NoTranspose", LU.rows, 1, &LU.elements[0], LU.rows, &pivots[0], &b[0], b.size(), &info);
assert(info == 0);
}
-void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b) {
- solveWithLUDecomposition(LU, pivots, &b[0], 1, b.size());
-}
-
// L (lower triangular) such that A = L L^T
// Inputs:
// - A : dim x dim symmetric matrix
@@ -282,60 +238,6 @@ void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stab
L.elements[i + j*dim] = 0;
}
-// 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 function dominates the running
-// time of this program.
-// 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 firstNonzeroIndex) {
- int dim = L.rows;
- Real c, s, x, y;
- for (int r = firstNonzeroIndex; r < dim; r++) {
- x = L.elt(r,r);
- y = *(v+r);
- Rrotg(&x, &y, &c, &s);
- Rrot(dim - r, &L.elements[r*(dim+1)], 1, v+r, 1, c, s);
- }
-}
-
-// Let lambdaGeometricMean be the geometric mean of the diagonal
-// elements of L. Whenever a diagonal element is less than
-// CHOLESKY_STABILIZE_THRESHOLD * lambdaGeometricMean, update the
-// cholesky decomposition L -> L' so that
-//
-// L' L'^T = L L^T + lambdaGeometricMean^2 u_i u_i^T
-//
-// where u is a unit vector in the i-th direction. The index i is
-// stored in the vector updateIndices.
-//
-void stabilizeCholesky(Matrix &L,
- Vector &updateVector,
- vector<int> &updateIndices,
- Real &lambdaGeometricMean) {
- int dim = L.rows;
- assert(L.cols == dim);
-
- double averageLogLambda = 0;
- for (int i = 0; i < dim; i++) {
- averageLogLambda += log(cast2double(L.elt(i,i)));
- }
- lambdaGeometricMean = Real(exp(averageLogLambda/dim));
-
- Real lambdaThreshold = CHOLESKY_STABILIZE_THRESHOLD * lambdaGeometricMean;
- for (int i = 0; i < dim; i++) {
- if (L.elt(i,i) < lambdaThreshold) {
- fillVector(updateVector, 0);
- updateVector[i] = lambdaGeometricMean;
- choleskyUpdate(L, &updateVector[0], i);
- updateIndices.push_back(i);
- }
- }
-}
-
void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb) {
int dim = L.rows;
assert(L.cols == dim);
@@ -376,45 +278,3 @@ void matrixSolveWithCholesky(Matrix &ACholesky, Matrix &X) {
lowerTriangularSolve(ACholesky, &X.elements[0], X.cols, X.rows);
lowerTriangularTransposeSolve(ACholesky, &X.elements[0], X.cols, X.rows);
}
-
-// Find the column with the maximum absolute value entry to the right
-// of (r,c) (inclusive)
-int maxColRightOf(Matrix &A, int r, int c) {
- int cMax = 0;
- Real aMax = 0;
- for (int k = c; k < A.cols; k++) {
- Real a = abs(A.elt(r,k));
- if (a > aMax) {
- cMax = k;
- aMax = a;
- }
- }
- return cMax;
-}
-
-// Return indices of a set of linearly independent rows of M, where
-// the basis matrix has eigenvalues whose absolute values exceed
-// thresh.
-vector<int> linearlyIndependentRowIndices(const Matrix &M) {
- vector<int> rows;
- Matrix A(M);
-
- int c = 0;
- for (int r = 0; r < A.rows && c < A.cols; r++) {
- A.swapCols(maxColRightOf(A, r, c), c);
- if (abs(A.elt(r, c)) > BASIC_ROW_THRESHOLD) {
- rows.push_back(r);
- // Add a multiple of column c to each column to the right so
- // that the entries to the right of (r,c) are zero.
- for (int k = c+1; k < A.cols; k++) {
- Real x = A.elt(r,k)/A.elt(r,c);
- for (int i = r; i < A.rows; i++)
- A.elt(i,k) -= x*A.elt(i,c);
- }
- c++;
- }
- }
- assert((int)rows.size() == A.cols);
-
- return rows;
-}
diff --git a/src/Matrix.h b/src/Matrix.h
index 5ffc80a..2b2cff8 100644
--- a/src/Matrix.h
+++ b/src/Matrix.h
@@ -121,15 +121,9 @@ void matrixScaleMultiplyAdd(Real alpha, Matrix &A, Matrix &B, Real beta, Matrix
// C := A*B
void matrixMultiply(Matrix &A, Matrix &B, Matrix &C);
-// B = A^T A
-void matrixSquare(Matrix &A, Matrix &B);
-
// Set block starting at (bRow, bCol) of B to A^T A
void matrixSquareIntoBlock(Matrix &A, Matrix &B, int bRow, int bCol);
-// A := L A L^T
-void lowerTriangularCongruence(Matrix &A, Matrix &L);
-
// A := L^{-1} A L^{-T}
void lowerTriangularInverseCongruence(Matrix &A, Matrix &L);
@@ -165,10 +159,6 @@ Real minEigenvalue(Matrix &A, Vector &workspace, Vector &eigenvalues);
// for use with 'solveWithLUDecomposition'
void LUDecomposition(Matrix &A, vector<Integer> &pivots);
-void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Real *B, int bcols, int ldb);
-
-void solveWithLUDecompositionTranspose(Matrix &LU, vector<Integer> &pivots, Real *B, int bcols, int ldb);
-
void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b);
// L (lower triangular) such that A = L L^T
@@ -179,27 +169,6 @@ void choleskyDecomposition(Matrix &A, Matrix &L);
void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas);
-// 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 function dominates the running
-// time of this program.
-// 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 firstNonzeroIndex);
-
-// Let lambdaGeometricMean be the geometric mean of the diagonal
-// elements of L. Whenever a diagonal element is less than
-// CHOLESKY_STABILIZE_THRESHOLD * lambdaGeometricMean, update the
-// cholesky decomposition L -> L' so that
-//
-// L' L'^T = L L^T + lambdaGeometricMean^2 u_i u_i^T
-//
-// where u is a unit vector in the i-th direction. The index i is
-// stored in the vector updateIndices.
-void stabilizeCholesky(Matrix &L, Vector &updateVector, vector<int> &updateIndices, Real &lambdaGeometricMean);
-
void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb);
void lowerTriangularSolve(Matrix &L, Vector &b);
@@ -214,9 +183,4 @@ void lowerTriangularTransposeSolve(Matrix &L, Vector &b);
// - X : dim x dim matrix
void matrixSolveWithCholesky(Matrix &ACholesky, Matrix &X);
-// Return indices of a set of linearly independent rows of M, where
-// the basis matrix has eigenvalues whose absolute values exceed
-// thresh.
-vector<int> linearlyIndependentRowIndices(const Matrix &M);
-
#endif // SDP_BOOTSTRAP_MATRIX_H_
diff --git a/src/SDPSolver.cpp b/src/SDPSolver.cpp
index 64cb3af..e1ca9ab 100644
--- a/src/SDPSolver.cpp
+++ b/src/SDPSolver.cpp
@@ -12,17 +12,14 @@ SDPSolver::SDPSolver(const SDP &sdp):
sdp(sdp),
x(sdp.primalObjective.size(), 0),
X(sdp.psdMatrixBlockDims()),
+ y(sdp.dualObjective.size(), 0),
Y(X),
dx(x),
dX(X),
+ dy(y),
dY(Y),
dualResidues(x),
- dualResiduesReduced(sdp.primalObjective.size() - sdp.dualObjective.size()),
PrimalResidues(X),
- FreeVarMatrixReduced(sdp.primalObjective.size() - sdp.dualObjective.size(), sdp.dualObjective.size()),
- FreeVarMatrixBasicLU(sdp.dualObjective.size(), sdp.dualObjective.size()),
- FreeVarMatrixBasicPivots(FreeVarMatrixBasicLU.rows),
- dualObjectiveReduced(sdp.dualObjective.size()),
XCholesky(X),
YCholesky(X),
Z(X),
@@ -46,51 +43,6 @@ SDPSolver::SDPSolver(const SDP &sdp):
eigenvaluesWorkspace.push_back(Vector(X.blocks[b].rows));
QRWorkspace.push_back(Vector(3*X.blocks[b].rows - 1));
}
-
- // Computations needed for free variable elimination
- basicIndices = linearlyIndependentRowIndices(sdp.FreeVarMatrix);
- for (int i = 0, p = 0; p < sdp.FreeVarMatrix.rows; p++)
- if (i < (int)basicIndices.size() && p == basicIndices[i])
- i++;
- else
- nonBasicIndices.push_back(p);
-
- // LU Decomposition of D_B
- for (int n = 0; n < FreeVarMatrixBasicLU.cols; n++)
- for (int m = 0; m < FreeVarMatrixBasicLU.rows; m++)
- FreeVarMatrixBasicLU.elt(m,n) = sdp.FreeVarMatrix.elt(basicIndices[m],n);
- LUDecomposition(FreeVarMatrixBasicLU, FreeVarMatrixBasicPivots);
-
- // Compute E = - D_N D_B^{-1}
- // ET = -D_N^T
- Matrix FreeVarMatrixReducedT(FreeVarMatrixReduced.cols, FreeVarMatrixReduced.rows);
- for (int p = 0; p < FreeVarMatrixReducedT.cols; p++)
- for (int n = 0; n < FreeVarMatrixReducedT.rows; n++)
- FreeVarMatrixReducedT.elt(n, p) = -sdp.FreeVarMatrix.elt(nonBasicIndices[p], n);
- // ET = D_B^{-1 T} ET = -D_B^{-1 T} D_N^T
- solveWithLUDecompositionTranspose(FreeVarMatrixBasicLU, FreeVarMatrixBasicPivots,
- &FreeVarMatrixReducedT.elements[0],
- FreeVarMatrixReducedT.cols,
- FreeVarMatrixReducedT.rows);
- // E = ET^T
- transpose(FreeVarMatrixReducedT, FreeVarMatrixReduced);
-
- // dualObjectiveReduced = D_B^{-T} f
- for (unsigned int n = 0; n < dualObjectiveReduced.size(); n++)
- dualObjectiveReduced[n] = sdp.dualObjective[n];
- solveWithLUDecompositionTranspose(FreeVarMatrixBasicLU,
- FreeVarMatrixBasicPivots,
- &dualObjectiveReduced[0], 1,
- dualObjectiveReduced.size());
-}
-
-Vector SDPSolver::freeVariableSolution() {
- Vector solution(sdp.dualObjective.size());
- for (unsigned int n = 0; n < solution.size(); n++) {
- solution[n] = dualResidues[basicIndices[n]];
- }
- solveWithLUDecomposition(FreeVarMatrixBasicLU, FreeVarMatrixBasicPivots, solution);
- return solution;
}
// result = b'^T a b', where b' = b \otimes 1
@@ -301,45 +253,8 @@ void computeSchurBlocks(const SDP &sdp,
}
}
-// x_B = g + E^T x_N
-void basicCompletion(const Vector &dualObjectiveReduced,
- const Matrix &FreeVarMatrixReduced,
- const vector<int> &basicIndices,
- const vector<int> &nonBasicIndices,
- Vector &x) {
- assert((int)basicIndices.size() == FreeVarMatrixReduced.cols);
- assert((int)nonBasicIndices.size() == FreeVarMatrixReduced.rows);
- assert((int)x.size() == FreeVarMatrixReduced.cols + FreeVarMatrixReduced.rows);
-
- #pragma omp parallel for schedule(static)
- for (unsigned int n = 0; n < basicIndices.size(); n++) {
- x[basicIndices[n]] = dualObjectiveReduced[n];
- for (unsigned int p = 0; p < nonBasicIndices.size(); p++)
- x[basicIndices[n]] += FreeVarMatrixReduced.elt(p, n) * x[nonBasicIndices[p]];
- }
-}
-
-// xReduced_N = x_N + E x_B
-void nonBasicShift(const Matrix &FreeVarMatrixReduced,
- const vector<int> &basicIndices,
- const vector<int> &nonBasicIndices,
- const Vector &x,
- Vector &xReduced) {
- assert((int)basicIndices.size() == FreeVarMatrixReduced.cols);
- assert((int)nonBasicIndices.size() == FreeVarMatrixReduced.rows);
- assert((int)x.size() == FreeVarMatrixReduced.cols + FreeVarMatrixReduced.rows);
- assert(nonBasicIndices.size() == xReduced.size());
-
- #pragma omp parallel for schedule(static)
- for (unsigned int p = 0; p < nonBasicIndices.size(); p++) {
- xReduced[p] = x[nonBasicIndices[p]];
- for (unsigned int n = 0; n < basicIndices.size(); n++)
- xReduced[p] += FreeVarMatrixReduced.elt(p,n) * x[basicIndices[n]];
- }
-}
-
void computeDualResidues(const SDP &sdp,
- const BlockDiagonalMatrix &Y,
+ const Vector &y,
const BlockDiagonalMatrix &BilinearPairingsY,
Vector &dualResidues) {
#pragma omp parallel for schedule(dynamic)
@@ -360,6 +275,9 @@ void computeDualResidues(const SDP &sdp,
dualResidues[p] -= BilinearPairingsY.blocks[*b].elt(ej_s+k, ej_r+k);
}
dualResidues[p] /= 2;
+
+ for (int n = 0; n < sdp.FreeVarMatrix.cols; n++)
+ dualResidues[p] -= sdp.FreeVarMatrix.elt(p, n)*y[n];
dualResidues[p] += sdp.primalObjective[p];
}
}
@@ -396,9 +314,11 @@ void constraintMatrixWeightedSum(const SDP &sdp, const Vector x, BlockDiagonalMa
}
void computeSchurRHS(const SDP &sdp,
- Vector &dualResidues,
- BlockDiagonalMatrix &Z,
- Vector &r) {
+ const Vector &dualResidues,
+ const BlockDiagonalMatrix &Z,
+ const Vector &x,
+ Vector &r,
+ Vector &s) {
for (unsigned int p = 0; p < r.size(); p++)
r[p] = -dualResidues[p];
@@ -418,6 +338,14 @@ void computeSchurRHS(const SDP &sdp,
}
}
}
+
+ #pragma omp parallel for schedule(static)
+ for (unsigned int n = 0; n < sdp.dualObjective.size(); n++) {
+ s[n] = sdp.dualObjective[n];
+ for (int p = 0; p < sdp.FreeVarMatrix.rows; p++) {
+ s[n] -= sdp.FreeVarMatrix.elt(p, n)*x[p];
+ }
+ }
}
// PrimalResidues = sum_p F_p x_p - X - F_0
@@ -430,18 +358,6 @@ void computePrimalResidues(const SDP &sdp,
PrimalResidues -= X;
}
-Real primalObjectiveValue(const SDP &sdp, const Vector &x) {
- return sdp.objectiveConst + dotProduct(sdp.primalObjective, x);
-}
-
-Real dualObjectiveValue(const SDP &sdp, const Vector &dualObjectiveReduced,
- const vector<int> &basicIndices, const Vector &dualResidues) {
- Real tmp = sdp.objectiveConst;
- for (unsigned int i = 0; i < dualObjectiveReduced.size(); i++)
- tmp += dualObjectiveReduced[i]*dualResidues[basicIndices[i]];
- return tmp;
-}
-
Real predictorCenteringParameter(const SDPSolverParameters ¶meters,
const bool isPrimalDualFeasible) {
return isPrimalDualFeasible ? Real(0) : parameters.infeasibleCenteringParameter;
@@ -577,13 +493,22 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
LUDecomposition(Q, Qpivots);
}
-void SDPSolver::solveSchurComplementEquation(Vector &dx) {
+
+// As inputs, dx and dy are the residues r_x and r_y on the right-hand
+// side of the Schur complement equation. As outputs, they are the
+// values for dx and dy.
+//
+void SDPSolver::solveSchurComplementEquation(Vector &dx, Vector &dy) {
// dx = SchurBlocksCholesky^{-1} dx
blockMatrixLowerTriangularSolve(SchurBlocksCholesky, dx);
basicKernelCoords.resize(Q.rows);
- // k_1 = -SchurUpdateLowRank^T dx
- vectorScaleMatrixMultiplyTransposeAdd(-1, SchurUpdateLowRank, dx, 0, basicKernelCoords);
+ // k_1 = r_y - SchurUpdateLowRank^T dx
+ for (unsigned int n = 0; n < dy.size(); n++)
+ basicKernelCoords[n] = dy[n];
+
+ vectorScaleMatrixMultiplyTransposeAdd(-1, SchurUpdateLowRank, dx, 1, basicKernelCoords);
+
// k_2 = -G^T dx
for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
int b = stabilizeBlockIndices[j];
@@ -615,6 +540,9 @@ void SDPSolver::solveSchurComplementEquation(Vector &dx) {
// dx = SchurBlocksCholesky^{-T} dx
blockMatrixLowerTriangularTransposeSolve(SchurBlocksCholesky, dx);
+ // dy = k_1
+ for (unsigned int n = 0; n < dy.size(); n++)
+ dy[n] = basicKernelCoords[n];
}
void SDPSolver::computeSearchDirection(const Real &beta,
@@ -632,11 +560,12 @@ void SDPSolver::computeSearchDirection(const Real &beta,
blockMatrixSolveWithCholesky(XCholesky, Z);
Z.symmetrize();
- // dx_k = -d_k + Tr(F_k Z)
- computeSchurRHS(sdp, dualResidues, Z, dx);
+ // (dx_RHS)_k = -d_k + Tr(F_k Z)
+ // dz_RHS = f - D^T x
+ computeSchurRHS(sdp, dualResidues, Z, x, dx, dy);
// dx_N = B_{NN}^{-1} dx_N, dx_B = E^T dx_N
- solveSchurComplementEquation(dx);
+ solveSchurComplementEquation(dx, dy);
// dX = R_p + sum_p F_p dx_p
constraintMatrixWeightedSum(sdp, dx, dX);
@@ -654,6 +583,7 @@ void SDPSolver::initialize(const SDPSolverParameters ¶meters) {
fillVector(x, 0);
X.setZero();
X.addDiagonal(parameters.initialMatrixScalePrimal);
+ fillVector(y, 0);
Y.setZero();
Y.addDiagonal(parameters.initialMatrixScaleDual);
}
@@ -667,26 +597,22 @@ SDPSolverTerminateReason SDPSolver::run(const SDPSolverParameters ¶meters,
if (timers["Run solver"].elapsed().wall >= parameters.maxRuntime * 1000000000LL)
return MaxRuntimeExceeded;
- // Maintain the invariant x_B = g + E^T x_N
- basicCompletion(dualObjectiveReduced, FreeVarMatrixReduced, basicIndices, nonBasicIndices, x);
-
choleskyDecomposition(X, XCholesky);
choleskyDecomposition(Y, YCholesky);
computeInvBilinearPairingsWithCholesky(XCholesky, sdp.bilinearBases, bilinearPairingsWorkspace, BilinearPairingsXInv);
computeBilinearPairings(Y, sdp.bilinearBases, bilinearPairingsWorkspace, BilinearPairingsY);
- // d_k = c_k - Tr(F_k Y)
- computeDualResidues(sdp, Y, BilinearPairingsY, dualResidues);
- nonBasicShift(FreeVarMatrixReduced, basicIndices, nonBasicIndices, dualResidues, dualResiduesReduced);
+ // d_k = c_k - Tr(F_k Y) - (D y)_k
+ computeDualResidues(sdp, y, BilinearPairingsY, dualResidues);
// PrimalResidues = sum_p F_p x_p - X - F_0 (F_0 is zero for now)
computePrimalResidues(sdp, x, X, PrimalResidues);
status.primalError = PrimalResidues.maxAbs();
- status.dualError = maxAbsVector(dualResiduesReduced);
- status.primalObjective = primalObjectiveValue(sdp, x);
- status.dualObjective = dualObjectiveValue(sdp, dualObjectiveReduced, basicIndices, dualResidues);
+ status.dualError = maxAbsVector(dualResidues);
+ status.primalObjective = sdp.objectiveConst + dotProduct(sdp.primalObjective, x);
+ status.dualObjective = sdp.objectiveConst + dotProduct(sdp.dualObjective, y);
const bool isPrimalFeasible = status.primalError < parameters.primalErrorThreshold;
const bool isDualFeasible = status.dualError < parameters.dualErrorThreshold;
@@ -694,8 +620,8 @@ SDPSolverTerminateReason SDPSolver::run(const SDPSolverParameters ¶meters,
if (isPrimalFeasible && isDualFeasible && isOptimal)
return PrimalDualOptimal;
- else if (isDualFeasible && status.dualObjective > parameters.maxDualObjective)
- return DualFeasibleMaxObjectiveExceeded;
+ else if (isDualFeasible && parameters.findDualFeasible)
+ return DualFeasible;
initializeSchurComplementSolver(BilinearPairingsXInv, BilinearPairingsY);
@@ -729,6 +655,8 @@ SDPSolverTerminateReason SDPSolver::run(const SDPSolverParameters ¶meters,
addVector(x, dx);
dX *= primalStepLength;
X += dX;
+ scaleVector(dy, dualStepLength);
+ addVector(y, dy);
dY *= dualStepLength;
Y += dY;
}
diff --git a/src/SDPSolver.h b/src/SDPSolver.h
index bf1bd7a..e8d3bd6 100644
--- a/src/SDPSolver.h
+++ b/src/SDPSolver.h
@@ -22,6 +22,7 @@ public:
int maxRuntime;
int checkpointInterval;
bool saveFinalCheckpoint;
+ bool findDualFeasible;
int precision;
int maxThreads;
Real dualityGapThreshold;
@@ -32,7 +33,6 @@ public:
Real feasibleCenteringParameter;
Real infeasibleCenteringParameter;
Real stepLengthReduction;
- Real maxDualObjective;
Real maxComplementarity;
void resetPrecision() {
@@ -44,7 +44,6 @@ public:
setPrecision(feasibleCenteringParameter, precision);
setPrecision(infeasibleCenteringParameter, precision);
setPrecision(stepLengthReduction, precision);
- setPrecision(maxDualObjective, precision);
setPrecision(maxComplementarity, precision);
}
@@ -53,7 +52,7 @@ public:
enum SDPSolverTerminateReason {
PrimalDualOptimal,
- DualFeasibleMaxObjectiveExceeded,
+ DualFeasible,
MaxComplementarityExceeded,
MaxIterationsExceeded,
MaxRuntimeExceeded,
@@ -84,26 +83,19 @@ public:
// current point
Vector x;
BlockDiagonalMatrix X;
+ Vector y;
BlockDiagonalMatrix Y;
// search direction
Vector dx;
BlockDiagonalMatrix dX;
+ Vector dy;
BlockDiagonalMatrix dY;
// discrepancies in dual and primal equality constraints
Vector dualResidues;
- Vector dualResiduesReduced;
BlockDiagonalMatrix PrimalResidues;
- // For free variable elimination
- Matrix FreeVarMatrixReduced;
- Matrix FreeVarMatrixBasicLU;
- vector<Integer> FreeVarMatrixBasicPivots;
- Vector dualObjectiveReduced;
- vector<int> basicIndices;
- vector<int> nonBasicIndices;
-
// intermediate computations
BlockDiagonalMatrix XCholesky;
BlockDiagonalMatrix YCholesky;
@@ -137,7 +129,7 @@ public:
SDPSolverTerminateReason run(const SDPSolverParameters ¶meters, const path checkpointFile);
void initializeSchurComplementSolver(const BlockDiagonalMatrix &BilinearPairingsXInv,
const BlockDiagonalMatrix &BilinearPairingsY);
- void solveSchurComplementEquation(Vector &dx);
+ void solveSchurComplementEquation(Vector &dx, Vector &dz);
void computeSearchDirection(const Real &beta, const Real &mu, const bool correctorPhase);
Vector freeVariableSolution();
void saveCheckpoint(const path &checkpointFile);
diff --git a/src/SDPSolverIO.cpp b/src/SDPSolverIO.cpp
index a974300..e36224c 100644
--- a/src/SDPSolverIO.cpp
+++ b/src/SDPSolverIO.cpp
@@ -52,6 +52,7 @@ ostream& operator<<(ostream& os, const SDPSolverParameters& p) {
os << "maxRuntime = " << p.maxRuntime << endl;
os << "checkpointInterval = " << p.checkpointInterval << endl;
os << "saveFinalCheckpoint = " << p.saveFinalCheckpoint << endl;
+ os << "findDualFeasible = " << p.findDualFeasible << endl;
os << "precision(actual) = " << p.precision << "(" << mpf_get_default_prec() << ")" << endl;
os << "maxThreads = " << p.maxThreads << endl;
os << "dualityGapThreshold = " << p.dualityGapThreshold << endl;
@@ -62,7 +63,6 @@ ostream& operator<<(ostream& os, const SDPSolverParameters& p) {
os << "feasibleCenteringParameter = " << p.feasibleCenteringParameter << endl;
os << "infeasibleCenteringParameter = " << p.infeasibleCenteringParameter << endl;
os << "stepLengthReduction = " << p.stepLengthReduction << endl;
- os << "maxDualObjective = " << p.maxDualObjective << endl;
os << "maxComplementarity = " << p.maxComplementarity << endl;
return os;
}
@@ -72,15 +72,15 @@ ostream &operator<<(ostream& os, const SDPSolverTerminateReason& r) {
case PrimalDualOptimal:
os << "found primal-dual optimal solution";
break;
+ case DualFeasible:
+ os << "found dual feasible solution";
+ break;
case MaxIterationsExceeded:
os << "maxIterations exceeded";
break;
case MaxRuntimeExceeded:
os << "maxRuntime exceeded";
break;
- case DualFeasibleMaxObjectiveExceeded:
- os << "found dual feasible solution with dualObjective exceeding maxDualObjective";
- break;
case MaxComplementarityExceeded:
os << "maxComplementarity exceeded";
break;
@@ -110,7 +110,7 @@ void SDPSolver::saveCheckpoint(const path &checkpointFile) {
boost::filesystem::ofstream ofs(checkpointFile);
boost::archive::text_oarchive ar(ofs);
cout << "Saving checkpoint to : " << checkpointFile << endl;
- boost::serialization::serializeSDPSolverState(ar, x, X, Y);
+ boost::serialization::serializeSDPSolverState(ar, x, X, y, Y);
timers["Save checkpoint"].start();
}
@@ -118,7 +118,7 @@ void SDPSolver::loadCheckpoint(const path &checkpointFile) {
boost::filesystem::ifstream ifs(checkpointFile);
boost::archive::text_iarchive ar(ifs);
cout << "Loading checkpoint from : " << checkpointFile << endl;
- boost::serialization::serializeSDPSolverState(ar, x, X, Y);
+ boost::serialization::serializeSDPSolverState(ar, x, X, y, Y);
}
void SDPSolver::saveSolution(const SDPSolverTerminateReason terminateReason, const path &outFile) {
@@ -132,8 +132,8 @@ void SDPSolver::saveSolution(const SDPSolverTerminateReason terminateReason, con
ofs << "primalError = " << status.primalError << ";\n";
ofs << "dualError = " << status.dualError << ";\n";
ofs << "runtime = " << float(timers["Run solver"].elapsed().wall)/1000000000 << ";\n";
- ofs << "freeVariables = " << freeVariableSolution() << ";\n";
+ ofs << "y = " << y << ";\n";
+ ofs << "Y = " << Y << ";\n";
ofs << "x = " << x << ";\n";
ofs << "X = " << X << ";\n";
- ofs << "Y = " << Y << ";\n";
}
diff --git a/src/main.cpp b/src/main.cpp
index 15eef50..81a8d1d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -117,6 +117,9 @@ int main(int argc, char** argv) {
("saveFinalCheckpoint",
po::value<bool>(¶meters.saveFinalCheckpoint)->default_value(true),
"Save a final checkpoint after terminating (useful to turn off when debugging).")
+ ("findDualFeasible",
+ po::value<bool>(¶meters.findDualFeasible)->default_value(false),
+ "Terminate once a dual feasible solution is found.")
("maxIterations",
po::value<int>(¶meters.maxIterations)->default_value(500),
"Maximum number of iterations to run the solver.")
@@ -155,9 +158,6 @@ int main(int argc, char** argv) {
po::value<Real>(¶meters.stepLengthReduction)->default_value(Real("0.7")),
"Shrink each newton step by this factor (smaller means slower, more "
"stable convergence). Corresponds to SDPA's gammaStar.")
- ("maxDualObjective",
- po::value<Real>(¶meters.maxDualObjective)->default_value(Real("1e10")),
- "Terminate if the dual objective exceeds this value.")
("maxComplementarity",
po::value<Real>(¶meters.maxComplementarity)->default_value(Real("1e100")),
"Terminate if the complementarity mu = Tr(X Y)/dim(X) exceeds this value.")
diff --git a/src/serialize.h b/src/serialize.h
index 376cacd..bfd24d2 100644
--- a/src/serialize.h
+++ b/src/serialize.h
@@ -57,9 +57,10 @@ namespace boost {
}
template<class Archive>
- void serializeSDPSolverState(Archive& ar, Vector &x, BlockDiagonalMatrix &X, BlockDiagonalMatrix &Y) {
+ void serializeSDPSolverState(Archive& ar, Vector &x, BlockDiagonalMatrix &X, Vector &y, BlockDiagonalMatrix &Y) {
ar & x;
ar & X;
+ ar & y;
ar & Y;
}
diff --git a/src/tests.cpp b/src/tests.cpp
index e99d96e..a0de85d 100644
--- a/src/tests.cpp
+++ b/src/tests.cpp
@@ -27,11 +27,6 @@ void testCholeskyStabilize() {
choleskyDecomposition(A,L);
cout << "A = " << A << ";\n";
cout << "L = " << L << ";\n";
- stabilizeCholesky(L, updateVector, updateIndices, lambdaGM);
- cout << "updateIndices = " << updateIndices << ";\n";
- cout << "L = " << L << "\n;";
- cout << "lambdaGM = " << lambdaGM << ";\n";
-
vector<Integer> stabilizeIndices;
vector<Real> stabilizeLambdas;
@@ -42,41 +37,6 @@ void testCholeskyStabilize() {
cout << "stabilizeLambdas = " << stabilizeLambdas << ";\n";
}
-void testLinearlyIndependentRowIndices() {
- Matrix A(8,3);
- A.elt(0,0)=0;
- A.elt(0,1)=0;
- A.elt(0,2)=0;
- A.elt(1,0)=1;
- A.elt(1,1)=0;
- A.elt(1,2)=1;
- A.elt(2,0)=0;
- A.elt(2,1)=1;
- A.elt(2,2)=0;
- A.elt(3,0)=0;
- A.elt(3,1)=0;
- A.elt(3,2)=0;
- A.elt(4,0)=0;
- A.elt(4,1)=1;
- A.elt(4,2)=1;
- A.elt(5,0)=1;
- A.elt(5,1)=1;
- A.elt(5,2)=0;
- A.elt(6,0)=1;
- A.elt(6,1)=0;
- A.elt(6,2)=1;
- A.elt(7,0)=1;
- A.elt(7,1)=0;
- A.elt(7,2)=1;
-
- cout << "A = " << A << ";\n";
-
- vector<int> rows = linearlyIndependentRowIndices(A);
-
- cout << "Aprime = " << A << ";\n";
- cout << "rows = " << rows << ";\n";
-}
-
void testCholeskyUpdate() {
Matrix A(4,4);
Matrix B(A);
--
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