[sdpb] 116/233: Commented SDP.h; changed some names
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 36a4f89eb390bcd742ad6699808a5407e1155218
Author: David Simmons-Duffin <davidsd at gmail.com>
Date: Fri Jan 2 19:13:09 2015 +1300
Commented SDP.h; changed some names
---
src/SDP.cpp | 78 +++++++++----------
src/SDP.h | 246 ++++++++++++++++++++++++++++++++++++++++++++++++++--------
src/parse.cpp | 51 ++++++------
3 files changed, 278 insertions(+), 97 deletions(-)
diff --git a/src/SDP.cpp b/src/SDP.cpp
index 6cc3f59..aa01dc8 100644
--- a/src/SDP.cpp
+++ b/src/SDP.cpp
@@ -24,69 +24,69 @@ Matrix sampleBilinearBasis(const int maxDegree,
return b;
}
-SampledMatrixPolynomial
-samplePolynomialVectorMatrix(const PolynomialVectorMatrix &m) {
- SampledMatrixPolynomial s;
+DualConstraintGroup
+polVecMatToDualConstraintGroup(const PolynomialVectorMatrix &m) {
+ DualConstraintGroup g;
assert(m.rows == m.cols);
- s.dim = m.rows;
- s.degree = m.degree();
+ g.dim = m.rows;
+ g.degree = m.degree();
- int numSamples = s.degree + 1;
- int numConstraints = numSamples * s.dim * (s.dim + 1)/2;
+ int numSamples = g.degree + 1;
+ int numConstraints = numSamples * g.dim * (g.dim + 1)/2;
int vectorDim = m.elt(0, 0).size();
// The first element of each vector multiplies the constant 1
- s.constraintConstants = Vector(numConstraints);
+ g.constraintConstants = Vector(numConstraints);
// The rest multiply decision variables
- s.constraintMatrix = Matrix(numConstraints, vectorDim - 1);
+ g.constraintMatrix = Matrix(numConstraints, vectorDim - 1);
int p = 0;
- for (int c = 0; c < s.dim; c++) {
- for (int r = c; r < s.dim; r++) {
+ for (int c = 0; c < g.dim; c++) {
+ for (int r = c; r < g.dim; r++) {
for (int k = 0; k < numSamples; k++) {
Real x = m.samplePoints[k];
Real scale = m.sampleScalings[k];
- s.constraintConstants[p] = scale*m.elt(r, c)[0](x);
+ g.constraintConstants[p] = scale*m.elt(r, c)[0](x);
for (int n = 1; n < vectorDim; n++)
- s.constraintMatrix.elt(p, n-1) = -scale*m.elt(r, c)[n](x);
+ g.constraintMatrix.elt(p, n-1) = -scale*m.elt(r, c)[n](x);
p++;
}
}
}
- int delta1 = s.degree/2;
- s.bilinearBases.push_back(sampleBilinearBasis(delta1, numSamples,
+ int delta1 = g.degree/2;
+ g.bilinearBases.push_back(sampleBilinearBasis(delta1, numSamples,
m.bilinearBasis,
m.samplePoints,
m.sampleScalings));
- int delta2 = (s.degree - 1)/2;
+ int delta2 = (g.degree - 1)/2;
if (delta2 >= 0)
- s.bilinearBases
+ g.bilinearBases
.push_back(sampleBilinearBasis(delta2, numSamples,
m.bilinearBasis,
m.samplePoints,
multiplyVectors(m.samplePoints,
m.sampleScalings)));
- return s;
+ return g;
}
-SDP bootstrapSDP(const Vector &objective,
+SDP sdpFromConstraintGroups(const Vector &objective,
const Real &objectiveConst,
- const vector<SampledMatrixPolynomial> &sampledMatrixPols) {
+ const vector<DualConstraintGroup> &dualConstraintGroups) {
SDP sdp;
sdp.dualObjective = objective;
sdp.objectiveConst = objectiveConst;
- for (vector<SampledMatrixPolynomial>::const_iterator s = sampledMatrixPols.begin();
- s != sampledMatrixPols.end();
- s++) {
- sdp.dimensions.push_back(s->dim);
- sdp.degrees.push_back(s->degree);
+ for (vector<DualConstraintGroup>::const_iterator g = dualConstraintGroups.begin();
+ g != dualConstraintGroups.end();
+ g++) {
+ sdp.dimensions.push_back(g->dim);
+ sdp.degrees.push_back(g->degree);
sdp.primalObjective.insert(sdp.primalObjective.end(),
s->constraintConstants.begin(),
s->constraintConstants.end());
@@ -94,22 +94,22 @@ SDP bootstrapSDP(const Vector &objective,
sdp.FreeVarMatrix = Matrix(sdp.primalObjective.size(), sdp.dualObjective.size());
int p = 0;
- for (vector<SampledMatrixPolynomial>::const_iterator s = sampledMatrixPols.begin();
- s != sampledMatrixPols.end();
- s++) {
+ for (vector<DualConstraintGroup>::const_iterator g = dualConstraintGroups.begin();
+ g != dualConstraintGroups.end();
+ g++) {
vector<int> blocks;
- for (vector<Matrix>::const_iterator b = s->bilinearBases.begin();
- b != s->bilinearBases.end();
+ for (vector<Matrix>::const_iterator b = g->bilinearBases.begin();
+ b != g->bilinearBases.end();
b++) {
- assert(b->cols == s->degree + 1);
+ assert(b->cols == g->degree + 1);
blocks.push_back(sdp.bilinearBases.size());
sdp.bilinearBases.push_back(*b);
}
sdp.blocks.push_back(blocks);
- for (int k = 0; k < s->constraintMatrix.rows; k++, p++)
- for (int n = 0; n < s->constraintMatrix.cols; n++)
- sdp.FreeVarMatrix.elt(p, n) = s->constraintMatrix.elt(k, n);
+ for (int k = 0; k < g->constraintMatrix.rows; k++, p++)
+ for (int n = 0; n < g->constraintMatrix.cols; n++)
+ sdp.FreeVarMatrix.elt(p, n) = g->constraintMatrix.elt(k, n);
}
assert(p == static_cast<int>(sdp.primalObjective.size()));
@@ -117,17 +117,17 @@ SDP bootstrapSDP(const Vector &objective,
return sdp;
}
-SDP bootstrapPolynomialSDP(const Vector &affineObjective,
- const vector<PolynomialVectorMatrix> &polVectorMatrices) {
- vector<SampledMatrixPolynomial> sampledMatrixPols;
+SDP bootstrapSDP(const Vector &affineObjective,
+ const vector<PolynomialVectorMatrix> &polVectorMatrices) {
+ vector<DualConstraintGroup> dualConstraintGroups;
for (vector<PolynomialVectorMatrix>::const_iterator m = polVectorMatrices.begin();
m != polVectorMatrices.end(); m++)
- sampledMatrixPols.push_back(samplePolynomialVectorMatrix(*m));
+ dualConstraintGroups.push_back(polVecMatToDualConstraintGroup(*m));
Vector objective = affineObjective;
objective.erase(objective.begin());
Real objectiveConst = affineObjective[0];
- return bootstrapSDP(objective, objectiveConst, sampledMatrixPols);
+ return sdpFromConstraintGroups(objective, objectiveConst, dualConstraintGroups);
}
diff --git a/src/SDP.h b/src/SDP.h
index 7fa2112..77995ba 100644
--- a/src/SDP.h
+++ b/src/SDP.h
@@ -22,32 +22,156 @@
using std::vector;
using std::ostream;
-// The constraint matrices A_* are labeled by tuples
-// (j_1,r_1,s_1,k_1), (j_2,r_2,s_2,k_2), ...
-// The corresponding IndexTuples are
-// IndexTuple(p=1,r=r_1,s=s_1,k=k_1), IndexTuple(p=2,r=r_2,s=s_2,k=k_2), ...
+// The class SDP encodes a semidefinite program of the following form
+//
+// Dual: maximize f + b.y over y,Y such that
+// Tr(A_p Y) + (B y)_p = c_p (0 <= p < P)
+// Y >= 0
+// Primal: minimize f + c.x over x,X such that
+// X = \sum_p A_p x_p - C
+// B^T x = b
+// X >= 0
+//
+// where the data of the SDP has the following structure
+//
+// - b,y are vectors of length N
+//
+// - c,x are vectors of length P
+//
+// - f is a constant that we add to the objective functions for
+// convenience. It has no effect on the running of our algorithm.
+//
+// - B is a P x N matrix, called the free variable matrix
+//
+// - X and Y are block diagonal:
+//
+// X = BlockDiagonal(X^(0), ..., X^(bMax-1))
+// Y = BlockDiagonal(Y^(0), ..., Y^(bMax-1))
+//
+// Let us define Block_b(M) as BlockDiagonal(0,...,0,M,0,...,0),
+// where M is in the b-th block. Then X and Y can be written
+//
+// X = \sum_{0<=b<bMax} Block_b(X^(b))
+// Y = \sum_{0<=b<bMax} Block_b(Y^(b))
+//
+// - The constraints labeled by 0 <= p < P are in 1-to-1
+// correspondence with tuples
+//
+// p <-> (j,r,s,k) where 0 <= j < J,
+// 0 <= s < m_j,
+// 0 <= r <= s,
+// 0 <= k <= d_j,
+//
+// We often interchange the index p with the tuple (j,r,s,k). The
+// constraint matrices A_p are given by
+//
+// A_(j,r,s,k) = \sum_{b \in blocks[j]}
+// Block_b(\chi_{b,k} \chi_{b,k}^T \otimes E^{rs}),
+//
+// where
+// - E^{rs} is the symmetrization of the m_j x m_j matrix with a
+// 1 in entry (r,s) and zeros elsewhere.
+// - \chi_{b,k} is a vector of length (delta_b+1)
+// - \otimes denotes a tensor product
+// - each block above thus has dimension (delta_b+1)*m_j.
+// - blocks[j_1] and blocks[j_2] are disjoint if j_1 != j_2. Thus,
+// the block index b determines a unique j, but not vice-versa.
+//
+
+// Instead of working with the 1 to 1 correspondence p <-> (j,r,s,k),
+// it is convenient to collect tuples with the same index j. This
+// gives a 1-to-many correspondence
+//
+// j <-> { list of IndexTuple(p,r,s,k) }
+//
+// IndexTuple is simply a named version of the 4-tuple (p,r,s,k). A
+// given IndexTuple uniquely specifies a constraint matrix A_p.
+//
class IndexTuple {
public:
- int p;
- int r;
- int s;
- int k;
+ int p; // overall index of the constraint
+ int r; // first index for E^{rs}
+ int s; // second index for E^{rs}
+ int k; // index for \chi_{b,k}
IndexTuple(int p, int r, int s, int k): p(p), r(r), s(s), k(k) {}
IndexTuple() {}
};
class SDP {
public:
+ // bilinearBases is a vector of Matrices encoding the \chi_{b,k}
+ // that enter the constraint matrices. Specifically, \chi_{b,k} are
+ // the columns of bilinearBases[b],
+ //
+ // bilinearBases[b].elt(m,k) = (\chi_{b,k})_m (0 <= b < bMax,
+ // 0 <= k <= d_j,
+ // 0 <= m <= delta_b)
+ //
vector<Matrix> bilinearBases;
+
+ // FreeVarMatrix = B, a PxN matrix
Matrix FreeVarMatrix;
+
+ // primalObjective = c, a vector of length P
Vector primalObjective;
+
+ // dualObjective = b, a vector of length N
Vector dualObjective;
+
+ // objectiveConst = f
Real objectiveConst;
+
+ // dimensions[j] = m_j (0 <= j < J)
vector<int> dimensions;
+
+ // degrees[j] = d_j (0 <= j < J)
vector<int> degrees;
+
+ // blocks gives the 1-to-many mapping
+ //
+ // blocks[j] = {b_{j1}, b_{j2}, ...} (0 <= j < J)
+ //
+ // entering the constraint matrices A_p. blocks[j1] and blocks[j2]
+ // are disjoint unless j1==j2.
+ //
vector<vector<int> > blocks;
+
+ // constraintIndices gives the 1-to-many mapping described above
+ //
+ // constraintIndices[j] = { IndexTuple(p,r,s,k)
+ // for 0 <= s < m_j,
+ // 0 <= r <= s,
+ // 0 <= k <= d_j,
+ // with p the overall constraint index }
+ //
+ // This allows us to loop through the constraints A_p associated to
+ // each j.
+ //
vector<vector<IndexTuple> > constraintIndices;
+ // create the mapping j -> { IndexTuple(p,r,s,k) } described above
+ //
+ void initializeConstraintIndices() {
+ int p = 0;
+ for (unsigned int j = 0; j < dimensions.size(); j++) {
+ constraintIndices.push_back(vector<IndexTuple>(0));
+
+ for (int s = 0; s < dimensions[j]; s++) {
+ for (int r = 0; r <= s; r++) {
+ for (int k = 0; k <= degrees[j]; k++) {
+ constraintIndices[j].push_back(IndexTuple(p, r, s, k));
+ p++;
+ }
+ }
+ }
+ }
+ assert(p == static_cast<int>(primalObjective.size()));
+ }
+
+ // Dimensions of the blocks of X,Y (0 <= b < bMax)
+ //
+ // psdMatrixBlockDims()[b] = (delta_b+1)*m_j = length(\chi_{b,*})*m_j
+ //
vector<int> psdMatrixBlockDims() const {
vector<int> dims;
for (unsigned int j = 0; j < dimensions.size(); j++)
@@ -57,6 +181,10 @@ class SDP {
return dims;
}
+ // Dimensions of the bilinear pairing matrices U^(b) and V^(b) (0 <= b < bMax)
+ //
+ // bilinearPairingBlockDims()[b] = (d_j + 1)*m_j
+ //
vector<int> bilinearPairingBlockDims() const {
vector<int> dims;
for (unsigned int j = 0; j < dimensions.size(); j++)
@@ -66,6 +194,11 @@ class SDP {
return dims;
}
+ // Dimensions of the blocks S^(j) of the Schur complement matrix:
+ //
+ // schurBlockDims()[j] = (d_j+1)*m_j*(m_j+1)/2
+ // = length(constraintIndices[j])
+ //
vector<int> schurBlockDims() const {
vector<int> dims;
for (unsigned int j = 0; j < dimensions.size(); j++)
@@ -73,23 +206,7 @@ class SDP {
return dims;
}
- void initializeConstraintIndices() {
- int p = 0;
- for (unsigned int j = 0; j < dimensions.size(); j++) {
- constraintIndices.push_back(vector<IndexTuple>(0));
-
- for (int s = 0; s < dimensions[j]; s++) {
- for (int r = 0; r <= s; r++) {
- for (int k = 0; k <= degrees[j]; k++) {
- constraintIndices[j].push_back(IndexTuple(p, r, s, k));
- p++;
- }
- }
- }
- }
- assert(p == static_cast<int>(primalObjective.size()));
- }
-
+ // Print an SDP, for debugging purposes
friend ostream& operator<<(ostream& os, const SDP& sdp) {
os << "SDP(bilinearBases = " << sdp.bilinearBases
<< ", FreeVarMatrix = " << sdp.FreeVarMatrix
@@ -103,22 +220,82 @@ class SDP {
}
};
-class SampledMatrixPolynomial {
+// DualConstraintGroup represents a set of constraints of the form
+//
+// Tr(A_p Y) + (B y)_p = c_p
+//
+// for a fixed j in the definition of SDP above. Here p corresponds to
+//
+// p <-> (r,s,k) where 0 <= s < dim,
+// 0 <= r <= s,
+// 0 <= k <= degree
+//
+// 0 <= p < (degree+1)*dim*(dim+1)/2 = P'
+//
+// The constraints of a full SDP can be thought of as a collection of
+// DualConstraintGroups labeled by 0<=j<J.
+//
+// DualConstraintGroup's are currently only used as an intermediate
+// data structure between the matrix polynomials defining an MPP and a
+// full SDP. By directly combining DualConstraintGroups into an SDP
+// using sdpFromConstraintGroups, it is possible to define slightly
+// more general optimization problems than those produced by
+// bootstrapSDP. Perhaps this level of generality will be useful in
+// the future.
+//
+class DualConstraintGroup {
public:
int dim;
int degree;
+ // constraintMatrix = B, a P'xN Matrix
Matrix constraintMatrix;
+ // constraintConstants = c, a vector of length P'
Vector constraintConstants;
+ // bilinearBases is a vector of Matrices encoding the \chi_{b,k}
+ // entering the constraint matrices A_p, as described
+ // above. `bilinearBases' here has the structure of
+ // `bilinearBases[j]' above for some fixed j.
vector<Matrix> bilinearBases;
};
+// Let M(x) be a matrix whose entries are vectors of polynomials:
+//
+// M(x) = ( \vec P^{00}(x) ... \vec P^{m0}(x) )
+// ( ... )
+// ( \vec P^{0n}(x) ... \vec P^{mn}(x) )
+//
+// where each vector has length N+1:
+//
+// \vec P^{rs}(x) = (P^{rs}_{-1}(x), P^{rs}_0, ... , P^{rs}_{N-1}(x))
+//
+// Consider a vector y = (y_0, ..., y_{N-1}) of length N, and let
+// (1,y) denote the vector of length N+1 whose components are 1, given
+// by the components of y. As explained in the manual, the constraint
+//
+// (1,y) . M(x) is positive semidefinite
+//
+// is equivalent to a DualConstraintGroup
+//
+// Tr(A_p Y) + (B y)_p = c_p
+//
+// A PolynomialVectorMatrix contains the data needed to construct this
+// DualConstraintGroup:
+//
class PolynomialVectorMatrix {
public:
- int rows;
- int cols;
+ int rows; // rows of M
+ int cols; // cols of M
+ // elements of M, in row-major order
vector<vector<Polynomial> > elements;
+
+ // A list of real numbers x_k (0 <= k <= degree(M)) at which to
+ // sample M(x) to construct the \chi_{b,k}.
vector<Real> samplePoints;
+ // A list of real numbers s_k (0 <= k <= degree(M)) to scale M(x_k)
+ // and the corresponding \chi_{b,k}.
vector<Real> sampleScalings;
+ // bilinearBasis[m] = q_m(x) (0 <= m <= degree/2), where q_m is a
+ // polynomial with degree deg(q_m) = m.
vector<Polynomial> bilinearBasis;
inline const vector<Polynomial>& elt(const int r, const int c) const {
@@ -129,6 +306,7 @@ class PolynomialVectorMatrix {
return elements[r + c*rows];
}
+ // The maximal degree of any of the components P^{rs}_n(x).
int degree() const {
int d = 0;
for (vector<vector<Polynomial> >::const_iterator e = elements.begin();
@@ -140,13 +318,13 @@ class PolynomialVectorMatrix {
}
};
-SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatrix &m);
+DualConstraintGroup polVecMatToDualConstraintGroup(const PolynomialVectorMatrix &m);
-SDP bootstrapSDP(const Vector &objective,
- const Real &objectiveConst,
- const vector<SampledMatrixPolynomial> &sampledMatrixPols);
+SDP sdpFromConstraintGroups(const Vector &objective,
+ const Real &objectiveConst,
+ const vector<DualConstraintGroup> &dualConstraintGroups);
-SDP bootstrapPolynomialSDP(const Vector &affineObjective,
- const vector<PolynomialVectorMatrix> &polVectorMatrices);
+SDP bootstrapSDP(const Vector &affineObjective,
+ const vector<PolynomialVectorMatrix> &polVectorMatrices);
#endif // SDPB_SDP_H_
diff --git a/src/parse.cpp b/src/parse.cpp
index 1c6586c..85b61b2 100644
--- a/src/parse.cpp
+++ b/src/parse.cpp
@@ -51,18 +51,6 @@ Matrix parseMatrix(XMLElement *xml) {
return m;
}
-SampledMatrixPolynomial parseSampledMatrixPolynomial(XMLElement *xml) {
- SampledMatrixPolynomial s;
- s.dim = parseInt(xml->FirstChildElement("dim"));
- s.degree = parseInt(xml->FirstChildElement("degree"));
- s.constraintMatrix = parseMatrix(xml->FirstChildElement("constraintMatrix"));
- s.constraintConstants = parseVector(xml->FirstChildElement("constraintConstants"));
- s.bilinearBases = parseMany("bilinearBasisMatrix",
- parseMatrix,
- xml->FirstChildElement("bilinearBases"));
- return s;
-}
-
Polynomial parsePolynomial(XMLElement *xml) {
Polynomial p;
p.coefficients = parseMany("coeff", parseReal, xml);
@@ -87,21 +75,36 @@ PolynomialVectorMatrix parsePolynomialVectorMatrix(XMLElement *xml) {
SDP parseBootstrapSDP(XMLElement *xml) {
return bootstrapSDP(parseVector(xml->FirstChildElement("objective")),
- parseReal(xml->FirstChildElement("objectiveConst")),
- parseMany("sampledMatrixPolynomial",
- parseSampledMatrixPolynomial,
- xml->FirstChildElement("sampledPositiveMatrices")));
-}
-
-SDP parseBootstrapPolynomialSDP(XMLElement *xml) {
- return bootstrapPolynomialSDP(parseVector(xml->FirstChildElement("objective")),
- parseMany("polynomialVectorMatrix",
- parsePolynomialVectorMatrix,
- xml->FirstChildElement("polynomialVectorMatrices")));
+ parseMany("polynomialVectorMatrix",
+ parsePolynomialVectorMatrix,
+ xml->FirstChildElement("polynomialVectorMatrices")));
}
SDP readBootstrapSDP(const path sdpFile) {
XMLDocument doc;
doc.LoadFile(sdpFile.c_str());
- return parseBootstrapPolynomialSDP(doc.FirstChildElement("sdp"));
+ return parseBootstrapSDP(doc.FirstChildElement("sdp"));
}
+
+// Reading SDPs defined by DualConstraintGroups is not necessary --
+// remove this functionality?
+
+// DualConstraintGroup parseDualConstraintGroup(XMLElement *xml) {
+// DualConstraintGroup g;
+// g.dim = parseInt(xml->FirstChildElement("dim"));
+// g.degree = parseInt(xml->FirstChildElement("degree"));
+// g.constraintMatrix = parseMatrix(xml->FirstChildElement("constraintMatrix"));
+// g.constraintConstants = parseVector(xml->FirstChildElement("constraintConstants"));
+// g.bilinearBases = parseMany("bilinearBasisMatrix",
+// parseMatrix,
+// xml->FirstChildElement("bilinearBases"));
+// return g;
+// }
+
+// SDP parseBootstrapSDP(XMLElement *xml) {
+// return sdpFromConstraintGroups(parseVector(xml->FirstChildElement("objective")),
+// parseReal(xml->FirstChildElement("objectiveConst")),
+// parseMany("dualConstraintGroup",
+// parseDualConstraintGroup,
+// xml->FirstChildElement("dualConstraintGroups")));
+// }
--
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