Skip to content

Commit

Permalink
Fixing bugs in the SparseMatrix implementations of the variant 2 Dant…
Browse files Browse the repository at this point in the history 
…zig selector and affine QP initialization (which effected both the affine QP and LP solvers)
  • Loading branch information
@poulson
poulson committed Dec 31, 2016
1 parent c3ddc44 commit 30575bd
Show file tree
Hide file tree
Showing 7 changed files with 115 additions and 40 deletions.
73 changes: 73 additions & 0 deletions examples/optimization/DS.cpp
View file
@@ -0,0 +1,73 @@
/*
Copyright (c) 2009-2016, Jack Poulson
All rights reserved.
This file is part of Elemental and is under the BSD 2-Clause License,
which can be found in the LICENSE file in the root directory, or at
http://opensource.org/licenses/BSD-2-Clause
*/
#include <El.hpp>

typedef double Real;

int
main( int argc, char* argv[] )
{
El::Environment env( argc, argv );

try
{
const El::Int m = El::Input("--m","height of matrix",100);
const El::Int n = El::Input("--n","width of matrix",200);
// TODO(poulson): Add options for controlling IPM
const El::Int maxIter =
El::Input("--maxIter","maximum # of iter's",500);
const Real lambda = El::Input("--lambda","DS parameter",0.5);
const bool display = El::Input("--display","display matrices?",false);
const bool print = El::Input("--print","print matrices",false);
El::ProcessInput();
El::PrintInputReport();

El::SparseMatrix<Real> A;
El::Matrix<Real> b, xTrue;
El::Identity( A, m, n );
El::Uniform( b, m, 1 );
if( print )
{
El::Print( A, "A" );
El::Print( b, "b" );
}
if( display )
El::Display( A, "A" );

El::lp::affine::Ctrl<Real> affineCtrl;
affineCtrl.mehrotraCtrl.print = true;

El::Matrix<Real> x;
El::Timer timer;
if( El::mpi::Rank() == 0 )
timer.Start();
El::DS( A, b, lambda, x, affineCtrl );
if( El::mpi::Rank() == 0 )
timer.Stop();
if( print )
El::Print( x, "x" );
const El::Int xZeroNorm = El::ZeroNorm( x );
if( El::mpi::Rank() == 0 )
{
El::Output("Dantzig Selector time: ",timer.Total()," secs");
El::Output("|| x ||_0 = ",xZeroNorm);
}
SoftThreshold( x, El::Sqrt(El::limits::Epsilon<Real>()) );
if( print )
El::Print( x, "xThresh" );
const El::Int xZeroNormThresh = El::ZeroNorm( x );
if( El::mpi::Rank() == 0 )
{
El::Output("|| xThresh ||_0 = ",xZeroNormThresh);
}
}
catch( std::exception& e ) { El::ReportException(e); }

return 0;
}
4 changes: 3 additions & 1 deletion examples/optimization/MPS.cpp
View file
Expand Up @@ -73,7 +73,9 @@ void SparseLoadAndSolve

timer.Start();
El::AffineLPSolution<El::Matrix<Real>> solution;
El::LP( problem, solution );
El::lp::affine::Ctrl<Real> ctrl;
ctrl.mehrotraCtrl.print = true;
El::LP( problem, solution, ctrl );
El::Output("Solving took ",timer.Stop()," seconds");
if( print )
{
Expand Down
3 changes: 2 additions & 1 deletion examples/optimization/RandomFeasibleLP.cpp
View file
Expand Up @@ -21,7 +21,8 @@ void RandomFeasibleLP( El::Int m, El::Int n, El::Int k )
El::Uniform( xFeas, n, 1 ); // Sample over B_1(0)
El::Uniform( sFeas, k, 1, Real(1), Real(1) ); // Sample over B_1(1)

El::Matrix<Real> A, G, b, c, h;
El::Matrix<Real> A, G;
El::Matrix<Real> b, c, h;
El::Uniform( A, m, n );
El::Uniform( G, k, n );
El::Gemv( El::NORMAL, Real(1), A, xFeas, b );
Expand Down
50 changes: 25 additions & 25 deletions src/optimization/models/DS.cpp
View file
Expand Up @@ -378,10 +378,34 @@ void Var2

// c := [1;1;0;0]
// ==============
Zeros( c, 3*n, 1 );
Zeros( c, 3*n+m, 1 );
auto cuv = c( IR(0,2*n), ALL );
Fill( cuv, Real(1) );

// G := | -I 0 0 0 |
// | 0 -I 0 0 |
// | 0 0 0 I |
// | 0 0 0 -I |
// ===================
Zeros( G, 4*n, 3*n+m );
G.Reserve( 4*n );
for( Int i=0; i<4*n; ++i )
{
if( i < 2*n )
G.QueueUpdate( i, i, Real(-1) );
else if( i < 3*n )
G.QueueUpdate( i, i+m, Real(+1) );
else
G.QueueUpdate( i, i+(m-n), Real(-1) );
}
G.ProcessQueues();

// h := [0;0;lambda e;lambda e]
// ============================
Zeros( h, 4*n, 1 );
auto ht = h( IR(2*n,4*n), ALL );
Fill( ht, lambda );

// \hat A := | A, -A, I, 0 |
// | 0, 0, A^T, I |
// ===========================
Expand All @@ -407,30 +431,6 @@ void Var2
auto b0 = bHat( IR(0,m), ALL );
b0 = b;

// G := | -I 0 0 0 |
// | 0 -I 0 0 |
// | 0 0 0 I |
// | 0 0 0 -I |
// ===================
Zeros( G, 4*n, 3*n+m );
G.Reserve( 4*n );
for( Int i=0; i<4*n; ++i )
{
if( i < 2*n )
G.QueueUpdate( i, i, Real(-1) );
else if( i < 3*n )
G.QueueUpdate( i, i+m, Real(+1) );
else
G.QueueUpdate( i, i+(m-n), Real(-1) );
}
G.ProcessQueues();

// h := [0;0;lambda e;lambda e]
// ============================
Zeros( h, 4*n, 1 );
auto ht = h( IR(2*n,4*n), ALL );
Fill( ht, lambda );

// Solve the affine LP
// ===================
Matrix<Real> xHat, y, z, s;
Expand Down
7 changes: 4 additions & 3 deletions src/optimization/solvers/LP/MPS.hpp
View file
Expand Up @@ -1017,7 +1017,7 @@ void Helper
if( compressed )
LogicError("Compressed MPS is not yet supported");
const bool upperBoundImplicitlyNonnegative = false;
const bool metadataSummary = false;
const bool metadataSummary = true;

MPSReader reader
( filename, compressed, minimize, upperBoundImplicitlyNonnegative );
Expand Down Expand Up @@ -1061,7 +1061,7 @@ void Helper
if( compressed )
LogicError("Compressed MPS is not yet supported");
const bool upperBoundImplicitlyNonnegative = false;
const bool metadataSummary = false;
const bool metadataSummary = true;

MPSReader reader
( filename, compressed, minimize, upperBoundImplicitlyNonnegative );
Expand Down Expand Up @@ -1102,14 +1102,15 @@ void Helper
if( compressed )
LogicError("Compressed MPS is not yet supported");
const bool upperBoundImplicitlyNonnegative = false;
const bool metadataSummary = false;
const bool metadataSummary = true;

MPSReader reader
( filename, compressed, minimize, upperBoundImplicitlyNonnegative );
const MPSMeta& meta = reader.Meta();
if( metadataSummary )
meta.PrintSummary();

Output("m=",meta.m,", n=",meta.n,", k=",meta.k);
Zeros( problem.c, meta.n, 1 );
Zeros( problem.A, meta.m, meta.n );
Zeros( problem.b, meta.m, 1 );
Expand Down
16 changes: 7 additions & 9 deletions src/optimization/solvers/LP/affine/IPM/Mehrotra.cpp
View file
Expand Up @@ -1083,10 +1083,6 @@ void EquilibratedMehrotra
{
try
{
J = JOrig;
J.FreezeSparsity();
UpdateDiagonal( J, Real(1), regTmp );

if( wMaxNorm >= ctrl.ruizEquilTol )
SymmetricRuizEquil( J, dInner, ctrl.ruizMaxIter, ctrl.print );
else if( wMaxNorm >= ctrl.diagEquilTol )
Expand Down Expand Up @@ -1257,6 +1253,9 @@ void EquilibratedMehrotra
residual.primalConic,
residual.dualConic,
solution.z, d );
J = JOrig;
J.FreezeSparsity();
UpdateDiagonal( J, Real(1), regTmp );

// Solve for the direction
// -----------------------
Expand Down Expand Up @@ -1501,11 +1500,6 @@ void EquilibratedMehrotra
{
try
{
J = JOrig;
J.FreezeSparsity();
J.LockedDistGraph().multMeta = JStatic.LockedDistGraph().multMeta;
UpdateDiagonal( J, Real(1), regTmp );

if( commRank == 0 && ctrl.time )
timer.Start();
if( wMaxNorm >= ctrl.ruizEquilTol )
Expand Down Expand Up @@ -1696,6 +1690,10 @@ void EquilibratedMehrotra
residual.primalConic,
residual.dualConic,
solution.z, d );
J = JOrig;
J.FreezeSparsity();
J.LockedDistGraph().multMeta = JStatic.LockedDistGraph().multMeta;
UpdateDiagonal( J, Real(1), regTmp );

// Solve for the direction
// -----------------------
Expand Down
2 changes: 1 addition & 1 deletion src/optimization/solvers/QP/affine/IPM/util/Initialize.cpp
View file
Expand Up @@ -329,7 +329,7 @@ void Initialize
{
EL_DEBUG_CSE
const Int m = b.Height();
const Int n = c.Width();
const Int n = c.Height();
const Int k = h.Height();
if( primalInit )
{
Expand Down

0 comments on commit 30575bd

Please sign in to comment.