Commit 392a0c6b authored by David Bommes's avatar David Bommes
Browse files

improvements of AQP and QuadraticProblem

parent bbb6d6da
Pipeline #2980 passed with stage
in 6 minutes and 6 seconds
......@@ -48,8 +48,16 @@ public:
/// Destructor
~AcceleratedQuadraticProxy() {}
// solve
int solve(NProblemInterface* _quadratic_problem, NProblemInterface* _nonlinear_problem, const SMatrixD& _A, const VectorD& _b)
// solve without linear constraints
int solve(NProblemInterface* _quadratic_problem, NProblemInterface* _nonlinear_problem, bool _update_factorization = true)
{
SMatrixD A(0,_quadratic_problem->n_unknowns());
VectorD b(VectorD::Index(0));
return solve(_quadratic_problem, _nonlinear_problem, A, b, _update_factorization);
}
// solve with linear constraints
int solve(NProblemInterface* _quadratic_problem, NProblemInterface* _nonlinear_problem, const SMatrixD& _A, const VectorD& _b, bool _update_factorization = true)
{
// time solution procedure
COMISO::StopWatch sw; sw.start();
......@@ -78,7 +86,8 @@ public:
const double theta = (1.0-std::sqrt(1.0/accelerate_))/(1.0+std::sqrt(1.0/accelerate_));
// pre-factorize linear system
pre_factorize(_quadratic_problem, _A, _b);
if(_update_factorization)
pre_factorize(_quadratic_problem, _A, _b);
int iter=0;
while( iter < max_iters_)
......@@ -206,7 +215,10 @@ protected:
if( fx_ls <= fx + alpha_ls_*t*gtdx )
{
_rel_df = 1.0-fx_ls/fx;
_rel_df = std::abs(1.0-fx_ls/fx);
std::cerr << "LS improved objective function " << fx << " -> " << fx_ls << std::endl;
return t;
}
else
......
......@@ -10,6 +10,10 @@
#include <Eigen/Sparse>
//== NAMESPACES ===============================================================
namespace COMISO {
// this problem optimizes the quadratic functional 0.5*x^T A x -x^t b + c
class QuadraticProblem : public COMISO::NProblemInterface
{
......@@ -18,14 +22,20 @@ public:
// Sparse Matrix Type
// typedef Eigen::DynamicSparseMatrix<double,Eigen::ColMajor> SMatrixNP;
QuadraticProblem(SMatrixNP& _A, Eigen::VectorXd _b, const double _c)
: A_(_A), c_(_c)
{
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
b_ = _b;
x_ = Eigen::VectorXd::Zero(A_.cols());
}
QuadraticProblem()
: A_(0,0), b_(Eigen::VectorXd::Index(0)), c_(0.0)
{
}
QuadraticProblem(SMatrixNP& _A, Eigen::VectorXd& _b, const double _c)
: A_(_A), c_(_c)
{
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
b_ = _b;
x_ = Eigen::VectorXd::Zero(A_.cols());
}
// number of unknowns
......@@ -77,6 +87,24 @@ public:
// advanced properties
virtual bool constant_hessian() { return true; }
void set_A(const SMatrixNP& _A)
{
A_ = _A;
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
x_ = Eigen::VectorXd::Zero(A_.cols());
}
void set_b(const Eigen::VectorXd& _b)
{
b_ = _b;
}
void set_c( const double _c)
{
c_ = _c;
}
private:
// quadratic problem 0.5*x^T A x -x^t b + c
......@@ -87,5 +115,8 @@ private:
Eigen::VectorXd x_;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // QUADRATICPROBLEM_H
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment