[sdpb] 83/233: Attempt at returning to polynomial SDPs -- cholesky woes on current test data
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:22 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 12ab4dacb5f412300903025c15dd917c01242694
Author: David Simmons-Duffin <dsd at athena.sns.ias.edu>
Date: Sat Sep 27 18:40:56 2014 -0400
Attempt at returning to polynomial SDPs -- cholesky woes on current test data
---
src/Polynomial.h | 66 ++++++++++++++++++++++++++++++++++++++++++++
src/SDP.cpp | 83 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
src/SDP.h | 12 ++++++++
src/Vector.h | 7 +++++
src/parse.cpp | 31 ++++++++++++++++++++-
5 files changed, 198 insertions(+), 1 deletion(-)
diff --git a/src/Polynomial.h b/src/Polynomial.h
new file mode 100644
index 0000000..112c70e
--- /dev/null
+++ b/src/Polynomial.h
@@ -0,0 +1,66 @@
+#ifndef SDP_BOOTSTRAP_POLYNOMIAL_H_
+#define SDP_BOOTSTRAP_POLYNOMIAL_H_
+
+#include "types.h"
+#include "util.h"
+#include "Vector.h"
+
+class Polynomial {
+public:
+ Vector coefficients;
+
+ Polynomial(): coefficients(1, 0) {}
+
+ int degree() const {
+ return coefficients.size() - 1;
+ };
+
+ Real operator()(const Real &x) const {
+ int deg = degree();
+ Real y = coefficients[deg];
+ for (int i = deg - 1; i >= 0; i--) {
+ y *= x;
+ y += coefficients[i];
+ }
+ return y;
+ }
+
+ friend ostream& operator<<(ostream& os, const Polynomial& p) {
+ for (int i = p.degree(); i >= 0; i--) {
+ os << p.coefficients[i];
+ if (i > 1)
+ os << "x^" << i << " + ";
+ else if (i == 1)
+ os << "x + ";
+ }
+ return os;
+ }
+
+};
+
+class PolynomialVectorMatrix {
+public:
+ int rows;
+ int cols;
+ vector<vector<Polynomial> > elements;
+
+ inline const vector<Polynomial>& elt(const int r, const int c) const {
+ return elements[r + c*rows];
+ }
+
+ inline vector<Polynomial>& elt(const int r, const int c) {
+ return elements[r + c*rows];
+ }
+
+ int degree() const {
+ int d = 0;
+ for (vector<vector<Polynomial> >::const_iterator e = elements.begin(); e != elements.end(); e++)
+ for (vector<Polynomial>::const_iterator p = e->begin(); p != e->end(); p++)
+ d = max(p->degree(), d);
+ return d;
+ }
+
+};
+
+
+#endif // SDP_BOOTSTRAP_POLYNOMIAL_H_
diff --git a/src/SDP.cpp b/src/SDP.cpp
index 73c12aa..4eb3174 100644
--- a/src/SDP.cpp
+++ b/src/SDP.cpp
@@ -1,5 +1,70 @@
#include "SDP.h"
+Matrix sampleBilinearBasis(const int basisSize,
+ const int numSamples,
+ const vector<Polynomial> &bilinearBasis,
+ const vector<Real> &samplePoints,
+ const vector<Real> &sampleScalings) {
+ Matrix b(basisSize, numSamples);
+ for (int k = 0; k < numSamples; k++) {
+ Real x = samplePoints[k];
+ Real scale = sqrt(sampleScalings[k]);
+ for (int i = 0; i < basisSize; i++)
+ b.elt(i, k) = scale*bilinearBasis[i](x);
+ }
+ return b;
+}
+
+SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatrix &m,
+ const vector<Polynomial> &bilinearBasis,
+ const vector<Real> &samplePoints,
+ const vector<Real> &sampleScalings) {
+ SampledMatrixPolynomial s;
+
+ assert(m.rows == m.cols);
+ s.dim = m.rows;
+ s.degree = m.degree();
+
+ int numSamples = s.degree + 1;
+ int numConstraints = numSamples * s.dim * (s.dim + 1)/2;
+ int vectorDim = m.elt(0,0).size();
+
+
+ // The first element of each vector multiplies the constant 1
+ s.constraintConstants = Vector(numConstraints);
+ // The rest multiply decision variables
+ s.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 k = 0; k < numSamples; k++) {
+ Real x = samplePoints[k];
+ Real scale = sampleScalings[k];
+
+ s.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);
+
+ p++;
+ }
+ }
+ }
+
+ int delta1 = s.degree/2;
+ s.bilinearBases.push_back(sampleBilinearBasis(delta1, numSamples, bilinearBasis,
+ samplePoints,
+ sampleScalings));
+ int delta2 = (s.degree - 1)/2;
+ if (delta2 >= 0)
+ s.bilinearBases.push_back(sampleBilinearBasis(delta2, numSamples, bilinearBasis,
+ samplePoints,
+ multiplyVectors(samplePoints, sampleScalings)));
+
+ return s;
+
+}
+
SDP bootstrapSDP(const Vector &objective,
const Real &objectiveConst,
const vector<SampledMatrixPolynomial> &sampledMatrixPols) {
@@ -42,3 +107,21 @@ SDP bootstrapSDP(const Vector &objective,
sdp.initializeConstraintIndices();
return sdp;
}
+
+SDP bootstrapPolynomialSDP(const Vector &affineObjective,
+ const vector<PolynomialVectorMatrix> &polVectorMatrices,
+ const vector<Polynomial> &bilinearBasis,
+ const vector<Real> &samplePoints,
+ const vector<Real> &sampleScalings) {
+ vector<SampledMatrixPolynomial> sampledMatrixPols;
+ for (vector<PolynomialVectorMatrix>::const_iterator m = polVectorMatrices.begin();
+ m != polVectorMatrices.end(); m++)
+ sampledMatrixPols.push_back(samplePolynomialVectorMatrix(*m, bilinearBasis, samplePoints, sampleScalings));
+
+ Vector objective = affineObjective;
+ objective.erase(objective.begin());
+ Real objectiveConst = affineObjective[0];
+
+ return bootstrapSDP(objective, objectiveConst, sampledMatrixPols);
+}
+
diff --git a/src/SDP.h b/src/SDP.h
index aa0ec13..7ca9957 100644
--- a/src/SDP.h
+++ b/src/SDP.h
@@ -8,6 +8,7 @@
#include "util.h"
#include "Vector.h"
#include "Matrix.h"
+#include "Polynomial.h"
using std::vector;
using std::ostream;
@@ -97,8 +98,19 @@ public:
vector<Matrix> bilinearBases;
};
+SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatrix &m,
+ const vector<Polynomial> &bilinearBasis,
+ const vector<Real> &samplePoints,
+ const vector<Real> &rescalings);
+
SDP bootstrapSDP(const Vector &objective,
const Real &objectiveConst,
const vector<SampledMatrixPolynomial> &sampledMatrixPols);
+SDP bootstrapPolynomialSDP(const Vector &affineObjective,
+ const vector<PolynomialVectorMatrix> &polVectorMatrices,
+ const vector<Polynomial> &bilinearBasis,
+ const vector<Real> &samplePoints,
+ const vector<Real> &sampleScalings);
+
#endif // SDP_BOOTSTRAP_SDP_H_
diff --git a/src/Vector.h b/src/Vector.h
index b27cf8e..04cfdcc 100644
--- a/src/Vector.h
+++ b/src/Vector.h
@@ -38,4 +38,11 @@ inline Real dotProduct(const Vector &u, const Vector v) {
return result;
}
+inline Vector multiplyVectors(const Vector &u, const Vector &v) {
+ Vector w(u.size());
+ for (unsigned int i = 0; i < w.size(); i++)
+ w[i] = u[i]*v[i];
+ return w;
+}
+
#endif // SDP_BOOTSTRAP_VECTOR_H_
diff --git a/src/parse.cpp b/src/parse.cpp
index 3caed01..088a0b4 100644
--- a/src/parse.cpp
+++ b/src/parse.cpp
@@ -2,6 +2,7 @@
#include "types.h"
#include "parse.h"
#include "tinyxml2.h"
+#include "Polynomial.h"
#include "SDP.h"
using std::vector;
@@ -51,6 +52,24 @@ SampledMatrixPolynomial parseSampledMatrixPolynomial(XMLElement *xml) {
return s;
}
+Polynomial parsePolynomial(XMLElement *polXml) {
+ Polynomial p;
+ p.coefficients = parseMany("coeff", parseReal, polXml);
+ return p;
+}
+
+vector<Polynomial> parsePolynomialVector(XMLElement *polVecXml) {
+ return parseMany("polynomial", parsePolynomial, polVecXml);
+}
+
+PolynomialVectorMatrix parsePolynomialVectorMatrix(XMLElement *polVecMatrixXml) {
+ PolynomialVectorMatrix m;
+ m.rows = parseInt(polVecMatrixXml->FirstChildElement("rows"));
+ m.cols = parseInt(polVecMatrixXml->FirstChildElement("cols"));
+ m.elements = parseMany("polynomialVector", parsePolynomialVector,
+ polVecMatrixXml->FirstChildElement("elements"));
+ return m;
+}
SDP parseBootstrapSDP(XMLElement *xml) {
return bootstrapSDP(parseVector(xml->FirstChildElement("objective")),
@@ -60,8 +79,18 @@ SDP parseBootstrapSDP(XMLElement *xml) {
xml->FirstChildElement("sampledPositiveMatrices")));
}
+SDP parseBootstrapPolynomialSDP(XMLElement *xml) {
+ return bootstrapPolynomialSDP(parseVector(xml->FirstChildElement("objective")),
+ parseMany("polynomialVectorMatrix",
+ parsePolynomialVectorMatrix,
+ xml->FirstChildElement("polynomialVectorMatrices")),
+ parsePolynomialVector(xml->FirstChildElement("bilinearBasis")),
+ parseVector(xml->FirstChildElement("samplePoints")),
+ parseVector(xml->FirstChildElement("sampleScalings")));
+}
+
SDP readBootstrapSDP(const path sdpFile) {
XMLDocument doc;
doc.LoadFile(sdpFile.c_str());
- return parseBootstrapSDP(doc.FirstChildElement("sdp"));
+ return parseBootstrapPolynomialSDP(doc.FirstChildElement("sdp"));
}
--
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