[sdpb] 118/233: Mostly commented SDPSolver.h
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:27 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 86a20010a0519c81de674487ef30e11e754caf3a
Author: David Simmons-Duffin <davidsd at gmail.com>
Date: Sat Jan 3 23:32:19 2015 +1300
Mostly commented SDPSolver.h
---
src/SDPSolver.cpp | 72 +++++++++++------------
src/SDPSolver.h | 172 ++++++++++++++++++++++++++++++++++++++++++++++++++----
2 files changed, 196 insertions(+), 48 deletions(-)
diff --git a/src/SDPSolver.cpp b/src/SDPSolver.cpp
index 326c0df..00c81d8 100644
--- a/src/SDPSolver.cpp
+++ b/src/SDPSolver.cpp
@@ -38,13 +38,13 @@ SDPSolver::SDPSolver(const SDP &sdp):
BilinearPairingsY(BilinearPairingsXInv),
SchurBlocks(sdp.schurBlockDims()),
SchurBlocksCholesky(SchurBlocks),
- SchurUpdateLowRank(sdp.FreeVarMatrix),
+ SchurOffDiagonal(sdp.FreeVarMatrix),
schurStabilizeIndices(SchurBlocks.blocks.size()),
schurStabilizeLambdas(SchurBlocks.blocks.size()),
stabilizeBlocks(SchurBlocks.blocks.size()),
Q(sdp.FreeVarMatrix.cols, sdp.FreeVarMatrix.cols),
Qpivots(sdp.FreeVarMatrix.cols),
- basicKernelCoords(Q.rows),
+ dyExtended(Q.rows),
StepMatrixWorkspace(X) {
// initialize bilinearPairingsWorkspace, eigenvaluesWorkspace, QRWorkspace
for (unsigned int b = 0; b < sdp.bilinearBases.size(); b++) {
@@ -426,32 +426,32 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
schurStabilizeLambdas,
cast2double(choleskyStabilizeThreshold));
timers["schurblocks/cholesky"].stop();
- // SchurUpdateLowRank = {{- 1, 0}, {E, G}}
+ // SchurOffDiagonal = {{- 1, 0}, {E, G}}
timers["make schur update"].resume();
- SchurUpdateLowRank.copyFrom(sdp.FreeVarMatrix);
- blockMatrixLowerTriangularSolve(SchurBlocksCholesky, SchurUpdateLowRank);
- int updateColumns = SchurUpdateLowRank.cols;
+ SchurOffDiagonal.copyFrom(sdp.FreeVarMatrix);
+ blockMatrixLowerTriangularSolve(SchurBlocksCholesky, SchurOffDiagonal);
+ int updateColumns = SchurOffDiagonal.cols;
stabilizeBlockIndices.resize(0);
stabilizeBlockUpdateRow.resize(0);
stabilizeBlockUpdateColumn.resize(0);
- for (unsigned int b = 0; b < SchurBlocks.blocks.size(); b++) {
- if (schurStabilizeIndices[b].size() > 0) {
- int startIndex = schurStabilizeIndices[b][0];
- int blockRows = SchurBlocks.blocks[b].rows - startIndex;
- int blockCols = schurStabilizeIndices[b].size();
+ for (unsigned int j = 0; j < SchurBlocks.blocks.size(); j++) {
+ if (schurStabilizeIndices[j].size() > 0) {
+ int startIndex = schurStabilizeIndices[j][0];
+ int blockRows = SchurBlocks.blocks[j].rows - startIndex;
+ int blockCols = schurStabilizeIndices[j].size();
- stabilizeBlockIndices.push_back(b);
- stabilizeBlockUpdateRow.push_back(SchurBlocks.blockStartIndices[b] + startIndex);
+ stabilizeBlockIndices.push_back(j);
+ stabilizeBlockUpdateRow.push_back(SchurBlocks.blockStartIndices[j] + startIndex);
stabilizeBlockUpdateColumn.push_back(updateColumns);
updateColumns += blockCols;
- stabilizeBlocks[b].setRowsCols(blockRows, blockCols);
- stabilizeBlocks[b].setZero();
- for (unsigned int c = 0; c < schurStabilizeIndices[b].size(); c++) {
- int r = schurStabilizeIndices[b][c] - startIndex;
- stabilizeBlocks[b].elt(r, c) = schurStabilizeLambdas[b][c];
+ stabilizeBlocks[j].setRowsCols(blockRows, blockCols);
+ stabilizeBlocks[j].setZero();
+ for (unsigned int c = 0; c < schurStabilizeIndices[j].size(); c++) {
+ int r = schurStabilizeIndices[j][c] - startIndex;
+ stabilizeBlocks[j].elt(r, c) = schurStabilizeLambdas[j][c];
}
}
}
@@ -470,11 +470,11 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
}
timers["make schur update"].stop();
timers["make Q"].resume();
- // Q = SchurUpdateLowRank^T SchurUpdateLowRank - {{0,0},{0,1}}
+ // Q = SchurOffDiagonal^T SchurOffDiagonal - {{0,0},{0,1}}
Q.setRowsCols(updateColumns, updateColumns);
Q.setZero();
- matrixSquareIntoBlock(SchurUpdateLowRank, Q, 0, 0);
+ matrixSquareIntoBlock(SchurOffDiagonal, Q, 0, 0);
// LowerRight(Q) = G^T G - 1
for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
@@ -493,13 +493,13 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
int r = stabilizeBlockUpdateColumn[j];
Rgemm("Transpose", "NoTranspose",
stabilizeBlocks[b].cols,
- SchurUpdateLowRank.cols,
+ SchurOffDiagonal.cols,
stabilizeBlocks[b].rows,
1,
&stabilizeBlocks[b].elements[0],
stabilizeBlocks[b].rows,
- &SchurUpdateLowRank.elt(p, 0),
- SchurUpdateLowRank.rows,
+ &SchurOffDiagonal.elt(p, 0),
+ SchurOffDiagonal.rows,
0,
&Q.elt(r, 0),
Q.rows);
@@ -507,8 +507,8 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
// UpperRight(Q) = LowerLeft(Q)^T
# pragma omp parallel for schedule(static)
- for (int c = 0; c < SchurUpdateLowRank.cols; c++)
- for (int r = SchurUpdateLowRank.cols; r < Q.rows; r++)
+ 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();
@@ -527,12 +527,12 @@ void SDPSolver::solveSchurComplementEquation(Vector &dx, Vector &dy) {
// dx = SchurBlocksCholesky^{-1} dx
blockMatrixLowerTriangularSolve(SchurBlocksCholesky, dx);
- basicKernelCoords.resize(Q.rows);
- // k_1 = r_y - SchurUpdateLowRank^T dx
+ dyExtended.resize(Q.rows);
+ // k_1 = r_y - SchurOffDiagonal^T dx
for (unsigned int n = 0; n < dy.size(); n++)
- basicKernelCoords[n] = dy[n];
+ dyExtended[n] = dy[n];
- vectorScaleMatrixMultiplyTransposeAdd(-1, SchurUpdateLowRank, dx, 1, basicKernelCoords);
+ vectorScaleMatrixMultiplyTransposeAdd(-1, SchurOffDiagonal, dx, 1, dyExtended);
// k_2 = -G^T dx
for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
@@ -541,17 +541,17 @@ void SDPSolver::solveSchurComplementEquation(Vector &dx, Vector &dy) {
int cTopLeft = stabilizeBlockUpdateColumn[j];
for (int c = 0; c < stabilizeBlocks[b].cols; c++) {
- basicKernelCoords[cTopLeft + c] = 0;
+ dyExtended[cTopLeft + c] = 0;
for (int r = 0; r < stabilizeBlocks[b].rows; r++)
- basicKernelCoords[cTopLeft + c] -= dx[pTopLeft + r] * stabilizeBlocks[b].elt(r, c);
+ dyExtended[cTopLeft + c] -= dx[pTopLeft + r] * stabilizeBlocks[b].elt(r, c);
}
}
// k = Q^{-1} k
- solveWithLUDecomposition(Q, Qpivots, basicKernelCoords);
+ solveWithLUDecomposition(Q, Qpivots, dyExtended);
- // dx = dx + SchurUpdateLowRank k_1
- vectorScaleMatrixMultiplyAdd(1, SchurUpdateLowRank, basicKernelCoords, 1, dx);
+ // dx = dx + SchurOffDiagonal k_1
+ vectorScaleMatrixMultiplyAdd(1, SchurOffDiagonal, dyExtended, 1, dx);
// dx += G k_2
for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
int b = stabilizeBlockIndices[j];
@@ -560,14 +560,14 @@ void SDPSolver::solveSchurComplementEquation(Vector &dx, Vector &dy) {
for (int c = 0; c < stabilizeBlocks[b].cols; c++)
for (int r = 0; r < stabilizeBlocks[b].rows; r++)
- dx[pTopLeft + r] += basicKernelCoords[cTopLeft + c] * stabilizeBlocks[b].elt(r, c);
+ dx[pTopLeft + r] += dyExtended[cTopLeft + c] * stabilizeBlocks[b].elt(r, c);
}
// dx = SchurBlocksCholesky^{-T} dx
blockMatrixLowerTriangularTransposeSolve(SchurBlocksCholesky, dx);
// dy = k_1
for (unsigned int n = 0; n < dy.size(); n++)
- dy[n] = basicKernelCoords[n];
+ dy[n] = dyExtended[n];
}
void SDPSolver::computeSearchDirection(const Real &beta,
diff --git a/src/SDPSolver.h b/src/SDPSolver.h
index 22941da..d145f44 100644
--- a/src/SDPSolver.h
+++ b/src/SDPSolver.h
@@ -24,6 +24,9 @@ using std::ostream;
using std::endl;
using boost::filesystem::path;
+// Parameters to control the behavior of an SDPSolver. See the manual
+// for a detailed description of each.
+//
class SDPSolverParameters {
public:
int maxIterations;
@@ -47,6 +50,10 @@ public:
Real choleskyStabilizeThreshold;
Real maxComplementarity;
+ // Set the precision of all Real parameters to equal 'precision'.
+ // This is necessary because 'precision' might be set (via the
+ // command line or a file) after initializing other parameters.
+ //
void resetPrecision() {
setPrecision(dualityGapThreshold, precision);
setPrecision(primalErrorThreshold, precision);
@@ -63,6 +70,9 @@ public:
friend ostream& operator<<(ostream& os, const SDPSolverParameters& p);
};
+// Reasons for terminating the solver. See the manual for a detailed
+// description of each.
+//
enum SDPSolverTerminateReason {
PrimalDualOptimal,
PrimalFeasible,
@@ -76,12 +86,15 @@ enum SDPSolverTerminateReason {
ostream &operator<<(ostream& os, const SDPSolverTerminateReason& r);
+// SDPSolverStatus contains the information needed to determine
+// whether to terminate or not.
+//
class SDPSolverStatus {
public:
- Real primalObjective;
- Real dualObjective;
- Real primalError;
- Real dualError;
+ Real primalObjective; // f + c . x
+ Real dualObjective; // f + b . y
+ Real primalError; // maxAbs(PrimalResidues)
+ Real dualError; // maxAbs(dualResidues)
Real dualityGap() const {
return abs(primalObjective - dualObjective) /
@@ -91,51 +104,187 @@ public:
friend ostream& operator<<(ostream& os, const SDPSolverStatus& s);
};
+// SDPSolver contains the data structures needed during the running of
+// the interior point algorithm. Each structure is allocated when an
+// SDPSolver is initialized, and reused in each iteration.
+//
class SDPSolver {
public:
+ // SDP to solve.
SDP sdp;
+
+ // Objective values and errors, re-evaluated each iteration
SDPSolverStatus status;
- // current point
+ /********************************************/
+ // The current point.
+
+ // a Vector of length P = sdp.primalObjective.size()
Vector x;
+
+ // a BlockDiagonalMatrix with structure given by sdp.psdMatrixBlockDims()
BlockDiagonalMatrix X;
+
+ // a Vector of length N = sdp.dualObjective.size()
Vector y;
+
+ // a BlockDiagonalMatrix with the same structure as X
BlockDiagonalMatrix Y;
- // search direction
+ /********************************************/
+ // The search direction.
+ //
+ // These quantities have the same structure as (x, X, y, Y). They
+ // are computed twice each iteration: once in the predictor step,
+ // and once in the corrector step.
+ //
Vector dx;
BlockDiagonalMatrix dX;
Vector dy;
BlockDiagonalMatrix dY;
- // discrepancies in dual and primal equality constraints
+ // Discrepancy in the dual equality constraints, a Vector of length
+ // P, called 'd' in the manual:
+ //
+ // dualResidues = c - Tr(A_* Y) - B y
+ //
Vector dualResidues;
+
+ // Discrepancy in the primal equality constraints, a
+ // BlockDiagonalMatrix with the same structure as X, called 'P' in
+ // the manual:
+ //
+ // PrimalResidues = \sum_P A_p x_p - X
+ //
BlockDiagonalMatrix PrimalResidues;
- // intermediate computations
+ /********************************************/
+ // Intermediate computations.
+
+ // the Cholesky decompositions of X and Y, each
+ // BlockDiagonalMatrices with the same structure as X and Y
BlockDiagonalMatrix XCholesky;
BlockDiagonalMatrix YCholesky;
+
+ // Z = X^{-1} (PrimalResidues Y - R), a BlockDiagonalMatrix with the
+ // same structure as X and Y
BlockDiagonalMatrix Z;
+
+ // R = mu I - X Y for the predictor step
+ // R = mu I - X Y - dX dY for the corrector step
+ // R is a BlockDiagonalMatrix with the same structure as X and Y.
BlockDiagonalMatrix R;
+
+ // Bilinear pairings needed for computing the Schur complement
+ // matrix. For example,
+ //
+ // BilinearPairingsXInv.blocks[b].elt(
+ // (d_j+1) s + k1,
+ // (d_j+1) r + k2
+ // ) = \chi_{b,k1}^T (X.blocks[b]^{-1})^{(s,r)} \chi_{b,k2}
+ //
+ // 0 <= k1,k2 <= sdp.degrees[j] = d_j
+ // 0 <= s,r < sdp.dimensions[j] = m_j
+ //
+ // where j corresponds to b and M^{(s,r)} denotes the (s,r)-th
+ // (d_j+1)x(d_j+1) block of M.
+ //
+ // BilinearPairingsXInv has one block for each block of X. The
+ // dimension of BilinearPairingsXInv.block[b] is (d_j+1)*m_j. See
+ // SDP.h for more information on d_j and m_j.
+ //
BlockDiagonalMatrix BilinearPairingsXInv;
+ //
+ // BilinearPairingsY is analogous to BilinearPairingsXInv, with
+ // X^{-1} -> Y.
+ //
BlockDiagonalMatrix BilinearPairingsY;
+
+ // The Schur complement matrix S: a BlockDiagonalMatrix with one
+ // block for each 0 <= j < J. SchurBlocks.blocks[j] has dimension
+ // (d_j+1)*m_j*(m_j+1)/2
+ //
BlockDiagonalMatrix SchurBlocks;
+
+ // SchurBlocksCholesky = L': the Cholesky decomposition of the
+ // stabilized Schur complement matrix S' = S + U U^T.
BlockDiagonalMatrix SchurBlocksCholesky;
- Matrix SchurUpdateLowRank;
+ // SchurOffDiagonal = L'^{-1} FreeVarMatrix: needed in solving the
+ // Schur complement equation.
+ Matrix SchurOffDiagonal;
+
+ // As explained in the manual, we use the `stabilized' Schur
+ // complement matrix S' = S + U U^T. Here, the 'update' matrix U
+ // has columns given by
+ //
+ // U = ( Lambda_{p_1} e_{p_1}, ..., Lambda_{p_M} e_{p_M} )
+ //
+ // where e_p is a unit vector in the p-th direction and the
+ // Lambda_{p_m} are constants. If p_m appears above, we say the
+ // direction p_m has been `stabilized.' Because U is sparse, we
+ // encode it in the smaller data structures below.
+ //
+ // Recall: S is block diagonal, with blocks labeled by 0 <= j < J.
+ //
+ // schurStabilizeIndices[j] = a list of which directions of
+ // SchurBlocks.blocks[j] have been stabilized (counting from 0 at
+ // the upper left of each block), for 0 <= j < J.
vector<vector<int> > schurStabilizeIndices;
+ //
+ // schurStabilizeLambdas[j] = a list of constants Lambda for each
+ // stabilized direction in schurStabilizeIndices[j], for 0 <= j < J.
vector<vector<Real> > schurStabilizeLambdas;
+ //
+ // a list of block indices {j_0, j_1, ..., j_{M-1} } for blocks
+ // which have at least one stabilized direction. We say the blocks
+ // j_m have been `stabilized.' Note that the number M of stabilized
+ // blocks could change with each iteration.
vector<int> stabilizeBlockIndices;
+ //
+ // For each 0 <= m < M, stabilizeBlockUpdateRow records the row of
+ // U corresponding to the first stabilized direction in the j_m'th
+ // block of SchurBlocks:
+ //
+ // stabilizeBlockUpdateRow[m] = schurStabilizeIndices[j_m][0] +
+ // SchurBlocks.blockStartIndices[j_m]
+ //
vector<int> stabilizeBlockUpdateRow;
+ //
+ // For each 0 <= m < M, stabilizeBlockUpdateColumn records the
+ // column of U corresponding to the first stabilized direction in
+ // the j_m'th block of SchurBlocks.
vector<int> stabilizeBlockUpdateColumn;
+ //
+ // For each 0 <= m < M, stabilizeBlocks is the submatrix of U
+ // corresponding to the j_m'th block of SchurBlocks.
vector<Matrix> stabilizeBlocks;
+ // Q = B' L'^{-T} L'^{-1} B' - {{0, 0}, {0, 1}}, where B' =
+ // (FreeVarMatrix U). Q is needed in the factorization of the Schur
+ // complement equation. Q has dimension N'xN', where
+ //
+ // N' = cols(B) + cols(U) = N + cols(U)
+ //
+ // where N is the dimension of the dual objective function. Note
+ // that N' could change with each iteration.
Matrix Q;
+ // a vector of length N', needed for the LU decomposition of Q.
vector<Integer> Qpivots;
- Vector basicKernelCoords;
- // additional workspace variables
+ // dyExtended = (dy z), where z are the extra coordinates introduced
+ // in the stabilized Schur complement equation.
+ Vector dyExtended;
+
+ /********************************************/
+ // Additional workspace variables
+
+ // A BlockDiagonalMatrix with the same structure as X, Y. Needed
+ // for computing step lengths which preserve positive
+ // semidefiniteness.
BlockDiagonalMatrix StepMatrixWorkspace;
+
+ //
vector<Matrix> bilinearPairingsWorkspace;
vector<Vector> eigenvaluesWorkspace;
vector<Vector> QRWorkspace;
@@ -148,7 +297,6 @@ public:
const Real &choleskyStabilizeThreshold);
void solveSchurComplementEquation(Vector &dx, Vector &dz);
void computeSearchDirection(const Real &beta, const Real &mu, const bool correctorPhase);
- Vector freeVariableSolution();
void saveCheckpoint(const path &checkpointFile);
void loadCheckpoint(const path &checkpointFile);
void saveSolution(const SDPSolverTerminateReason, const path &outFile);
--
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