[sdpb] 122/233: Made dualityGap a variable like primal/dual objectives
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:28 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 68e84ebf96eccf9902a4c5600cf70546bf962cfa
Author: David Simmons-Duffin <davidsd at gmail.com>
Date: Fri Jan 9 13:01:40 2015 +1300
Made dualityGap a variable like primal/dual objectives
---
src/SDPSolver.cpp | 43 ++++++++++++++++++-------------------------
src/SDPSolver.h | 11 +++++------
src/SDPSolverIO.cpp | 5 +++--
src/main.cpp | 2 +-
4 files changed, 27 insertions(+), 34 deletions(-)
diff --git a/src/SDPSolver.cpp b/src/SDPSolver.cpp
index 88e2b4f..4a4b7f7 100644
--- a/src/SDPSolver.cpp
+++ b/src/SDPSolver.cpp
@@ -46,7 +46,8 @@ SDPSolver::SDPSolver(const SDP &sdp, const SDPSolverParameters ¶meters):
Q(sdp.FreeVarMatrix.cols, sdp.FreeVarMatrix.cols),
Qpivots(sdp.FreeVarMatrix.cols),
dyExtended(Q.rows),
- StepMatrixWorkspace(X) {
+ StepMatrixWorkspace(X)
+{
// initialize bilinearPairingsWorkspace, eigenvaluesWorkspace, QRWorkspace
for (unsigned int b = 0; b < sdp.bilinearBases.size(); b++) {
bilinearPairingsWorkspace.push_back(Matrix(X.blocks[b].rows,
@@ -425,15 +426,13 @@ Real stepLength(BlockDiagonalMatrix &XCholesky,
void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &BilinearPairingsXInv,
const BlockDiagonalMatrix &BilinearPairingsY,
const Real &choleskyStabilizeThreshold) {
- timers["schurblocks/cholesky"].resume();
computeSchurBlocks(sdp, BilinearPairingsXInv, BilinearPairingsY, SchurBlocks);
choleskyDecompositionStabilized(SchurBlocks, SchurBlocksCholesky,
schurStabilizeIndices,
schurStabilizeLambdas,
cast2double(choleskyStabilizeThreshold));
- timers["schurblocks/cholesky"].stop();
+
// SchurOffDiagonal = {{- 1, 0}, {E, G}}
- timers["make schur update"].resume();
SchurOffDiagonal.copyFrom(sdp.FreeVarMatrix);
blockMatrixLowerTriangularSolve(SchurBlocksCholesky, SchurOffDiagonal);
int updateColumns = SchurOffDiagonal.cols;
@@ -474,8 +473,7 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
&stabilizeBlocks[b].elt(0, 0),
stabilizeBlocks[b].rows);
}
- timers["make schur update"].stop();
- timers["make Q"].resume();
+
// Q = SchurOffDiagonal^T SchurOffDiagonal - {{0,0},{0,1}}
Q.setRowsCols(updateColumns, updateColumns);
Q.setZero();
@@ -516,12 +514,9 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
for (int c = 0; c < SchurOffDiagonal.cols; c++)
for (int r = SchurOffDiagonal.cols; r < Q.rows; r++)
Q.elt(c, r) = Q.elt(r, c);
- timers["make Q"].stop();
- timers["LU of Q"].resume();
Qpivots.resize(Q.rows);
LUDecomposition(Q, Qpivots);
- timers["LU of Q"].stop();
}
@@ -623,19 +618,19 @@ SDPSolverTerminateReason SDPSolver::run(const path checkpointFile) {
primalObjective = sdp.objectiveConst + dotProduct(sdp.primalObjective, x);
dualObjective = sdp.objectiveConst + dotProduct(sdp.dualObjective, y);
+ dualityGap = abs(primalObjective - dualObjective) /
+ max(Real(abs(primalObjective) + abs(dualObjective)), Real(1));
- timers["cholesky"].resume();
choleskyDecomposition(X, XCholesky);
choleskyDecomposition(Y, YCholesky);
- timers["cholesky"].stop();
- timers["bilinear pairings"].resume();
+
computeInvBilinearPairingsWithCholesky(XCholesky, sdp.bilinearBases,
bilinearPairingsWorkspace,
BilinearPairingsXInv);
computeBilinearPairings(Y, sdp.bilinearBases,
bilinearPairingsWorkspace,
BilinearPairingsY);
- timers["bilinear pairings"].stop();
+
// d_k = c_k - Tr(F_k Y) - (D y)_k
computeDualResidues(sdp, y, BilinearPairingsY, dualResidues);
dualError = maxAbsVector(dualResidues);
@@ -644,9 +639,9 @@ SDPSolverTerminateReason SDPSolver::run(const path checkpointFile) {
computePrimalResidues(sdp, x, X, PrimalResidues);
primalError = PrimalResidues.maxAbs();
- const bool isPrimalFeasible = primalError < parameters.primalErrorThreshold;
- const bool isDualFeasible = dualError < parameters.dualErrorThreshold;
- const bool isOptimal = dualityGap() < parameters.dualityGapThreshold;
+ const bool isPrimalFeasible = primalError < parameters.primalErrorThreshold;
+ const bool isDualFeasible = dualError < parameters.dualErrorThreshold;
+ const bool isOptimal = dualityGap < parameters.dualityGapThreshold;
if (isPrimalFeasible && isDualFeasible && isOptimal)
return PrimalDualOptimal;
@@ -662,39 +657,37 @@ SDPSolverTerminateReason SDPSolver::run(const path checkpointFile) {
// functions for the current point
else if (iteration > parameters.maxIterations)
return MaxIterationsExceeded;
- timers["initialize schur solver"].resume();
+
initializeSchurComplementSolver(BilinearPairingsXInv, BilinearPairingsY,
parameters.choleskyStabilizeThreshold);
- timers["initialize schur solver"].stop();
Real mu = frobeniusProductSymmetric(X, Y)/X.dim;
if (mu > parameters.maxComplementarity)
return MaxComplementarityExceeded;
- timers["predictor"].resume();
+
// Mehrotra predictor solution for (dx, dX, dY)
Real betaPredictor = predictorCenteringParameter(parameters,
isPrimalFeasible && isDualFeasible);
computeSearchDirection(betaPredictor, mu, false);
- timers["predictor"].stop();
- timers["corrector"].resume();
+
// Mehrotra corrector solution for (dx, dX, dY)
Real betaCorrector = correctorCenteringParameter(parameters, X, dX, Y, dY, mu,
isPrimalFeasible && isDualFeasible);
computeSearchDirection(betaCorrector, mu, true);
- timers["corrector"].stop();
- timers["step length"].resume();
+
// Step length to preserve positive definiteness
primalStepLength = stepLength(XCholesky, dX, StepMatrixWorkspace,
eigenvaluesWorkspace, QRWorkspace, parameters);
dualStepLength = stepLength(YCholesky, dY, StepMatrixWorkspace,
eigenvaluesWorkspace, QRWorkspace, parameters);
- timers["step length"].stop();
+
if (isPrimalFeasible && isDualFeasible) {
primalStepLength = min(primalStepLength, dualStepLength);
dualStepLength = primalStepLength;
}
- printSolverInfo(iteration, mu, primalObjective, dualObjective, primalError, dualError, primalStepLength, dualStepLength,
+ printSolverInfo(iteration, mu, primalObjective, dualObjective, dualityGap,
+ primalError, dualError, primalStepLength, dualStepLength,
betaCorrector, sdp.dualObjective.size(), Q.rows);
// Update current point
diff --git a/src/SDPSolver.h b/src/SDPSolver.h
index 3c9aef2..0341b07 100644
--- a/src/SDPSolver.h
+++ b/src/SDPSolver.h
@@ -132,6 +132,7 @@ public:
//
Real primalObjective; // f + c . x
Real dualObjective; // f + b . y
+ Real dualityGap; // normalized difference of objectives
// Discrepancy in the primal equality constraints, a
// BlockDiagonalMatrix with the same structure as X, called 'P' in
@@ -290,15 +291,12 @@ public:
/********************************************/
// Methods
- // Create a new solver for a given SDP
+ // Create a new solver for a given SDP, with the given parameters
SDPSolver(const SDP &sdp, const SDPSolverParameters ¶meters);
- Real dualityGap() const {
- return abs(primalObjective - dualObjective) /
- max(Real(abs(primalObjective) + abs(dualObjective)), Real(1));
- }
-
+ // Run the solver, backing up to checkpointFile
SDPSolverTerminateReason run(const path checkpointFile);
+
void initializeSchurComplementSolver(const BlockDiagonalMatrix &BilinearPairingsXInv,
const BlockDiagonalMatrix &BilinearPairingsY,
const Real &choleskyStabilizeThreshold);
@@ -314,6 +312,7 @@ void printSolverInfo(int iteration,
Real mu,
Real primalObjective,
Real dualObjective,
+ Real dualityGap,
Real primalError,
Real dualError,
Real primalStepLength,
diff --git a/src/SDPSolverIO.cpp b/src/SDPSolverIO.cpp
index 5b4c603..6e85093 100644
--- a/src/SDPSolverIO.cpp
+++ b/src/SDPSolverIO.cpp
@@ -33,6 +33,7 @@ void printSolverInfo(int iteration,
Real mu,
Real primalObjective,
Real dualObjective,
+ Real dualityGap,
Real primalError,
Real dualError,
Real primalStepLength,
@@ -51,7 +52,7 @@ void printSolverInfo(int iteration,
mu.get_mpf_t(),
primalObjective.get_mpf_t(),
dualObjective.get_mpf_t(),
- dualityGap().get_mpf_t(),
+ dualityGap.get_mpf_t(),
primalError.get_mpf_t(),
dualError.get_mpf_t(),
primalStepLength.get_mpf_t(),
@@ -149,7 +150,7 @@ void SDPSolver::saveSolution(const SDPSolverTerminateReason terminateReason, con
ofs << "terminateReason = \"" << terminateReason << "\";\n";
ofs << "primalObjective = " << primalObjective << ";\n";
ofs << "dualObjective = " << dualObjective << ";\n";
- ofs << "dualityGap = " << dualityGap() << ";\n";
+ ofs << "dualityGap = " << dualityGap << ";\n";
ofs << "primalError = " << primalError << ";\n";
ofs << "dualError = " << dualError << ";\n";
ofs << "runtime = " << runtime << ";\n";
diff --git a/src/main.cpp b/src/main.cpp
index 81fde43..ca3698f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -68,7 +68,7 @@ int solveSDP(const path &sdpFile,
cout << endl;
cout << "primalObjective = " << solver.primalObjective << endl;
cout << "dualObjective = " << solver.dualObjective << endl;
- cout << "dualityGap = " << solver.dualityGap() << endl;
+ cout << "dualityGap = " << solver.dualityGap << endl;
cout << "primalError = " << solver.primalError << endl;
cout << "dualError = " << solver.dualError << endl;
cout << endl;
--
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