Commit baa6cb92 authored by Tobias Lasser's avatar Tobias Lasser

OptimizationProblem is now the new base class 'Problem' for problems, added...

OptimizationProblem is now the new base class 'Problem' for problems, added QuadricProblem incl. conversions between different types of problems, added CG solver
parent be7c5e6b
Pipeline #170242 passed with stages
in 22 minutes and 38 seconds
......@@ -9,17 +9,18 @@ Problem
.. doxygenclass:: elsa::Problem
RegularizationTerm
==================
.. doxygenclass:: elsa::RegularizationTerm
WLSProblem
==========
.. doxygenclass:: elsa::WLSProblem
OptimizationProblem
QuadricProblem
===================
.. doxygenclass:: elsa::OptimizationProblem
.. doxygenclass:: elsa::QuadricProblem
RegularizationTerm
==================
.. doxygenclass:: elsa::RegularizationTerm
......@@ -9,6 +9,11 @@ Solver
.. doxygenclass:: elsa::Solver
CG
==
.. doxygenclass:: elsa::CG
GradientDescent
===============
......
......@@ -62,7 +62,7 @@ namespace elsa
template <typename data_t>
DataContainer<data_t> LinearOperator<data_t>::apply(const DataContainer<data_t>& x)
DataContainer<data_t> LinearOperator<data_t>::apply(const DataContainer<data_t>& x) const
{
DataContainer<data_t> result(*_rangeDescriptor);
apply(x, result);
......@@ -70,13 +70,13 @@ namespace elsa
}
template <typename data_t>
void LinearOperator<data_t>::apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax)
void LinearOperator<data_t>::apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const
{
_apply(x, Ax);
}
template <typename data_t>
void LinearOperator<data_t>::_apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax)
void LinearOperator<data_t>::_apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const
{
if (_isLeaf) {
if (_isAdjoint) {
......@@ -132,7 +132,7 @@ namespace elsa
template <typename data_t>
DataContainer<data_t> LinearOperator<data_t>::applyAdjoint(const DataContainer<data_t>& y)
DataContainer<data_t> LinearOperator<data_t>::applyAdjoint(const DataContainer<data_t>& y) const
{
DataContainer<data_t> result(*_domainDescriptor);
applyAdjoint(y, result);
......@@ -140,13 +140,13 @@ namespace elsa
}
template <typename data_t>
void LinearOperator<data_t>::applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty)
void LinearOperator<data_t>::applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) const
{
_applyAdjoint(y, Aty);
}
template <typename data_t>
void LinearOperator<data_t>::_applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty)
void LinearOperator<data_t>::_applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) const
{
if (_isLeaf) {
if (_isAdjoint) {
......
......@@ -66,7 +66,7 @@ namespace elsa
*
* Please note: this method uses apply(x, Ax) to perform the actual operation.
*/
DataContainer<data_t> apply(const DataContainer<data_t>& x);
DataContainer<data_t> apply(const DataContainer<data_t>& x) const;
/**
* \brief apply the operator A to an element in the operator's domain
......@@ -78,7 +78,7 @@ namespace elsa
* classes. (Why is this method not virtual itself? Because you cannot have a non-virtual
* function overloading a virtual one [apply with one vs. two arguments]).
*/
void apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax);
void apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const;
/**
* \brief apply the adjoint of operator A to an element of the operator's range
......@@ -90,7 +90,7 @@ namespace elsa
*
* Please note: this method uses applyAdjoint(y, Aty) to perform the actual operation.
*/
DataContainer<data_t> applyAdjoint(const DataContainer<data_t>& y);
DataContainer<data_t> applyAdjoint(const DataContainer<data_t>& y) const;
/**
* \brief apply the adjoint of operator A to an element of the operator's range
......@@ -102,7 +102,7 @@ namespace elsa
* derived classes. (Why is this method not virtual itself? Because you cannot have a
* non-virtual function overloading a virtual one [applyAdjoint with one vs. two args]).
*/
void applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty);
void applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) const;
/// friend operator+ to support composition of LinearOperators (and its derivatives)
......@@ -141,10 +141,10 @@ namespace elsa
/// the apply method that has to be overridden in derived classes
virtual void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax);
virtual void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const;
/// the applyAdjoint method that has to be overridden in derived classes
virtual void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty);
virtual void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) const;
private:
/// pointers to nodes in the evaluation tree
......
......@@ -21,11 +21,11 @@ public:
: LinearOperator<data_t>(domain, range) {}
protected:
void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) override {
void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override {
Ax = 1;
}
void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) override {
void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) const override {
Aty = 3;
}
......
......@@ -46,7 +46,7 @@ namespace elsa
}
template <typename data_t>
LinearOperator<data_t>& LinearResidual<data_t>::getOperator() const
const LinearOperator<data_t>& LinearResidual<data_t>::getOperator() const
{
if (!_operator)
throw std::runtime_error("LinearResidual::getOperator: operator not present");
......@@ -55,7 +55,7 @@ namespace elsa
}
template <typename data_t>
const DataContainer<data_t>& LinearResidual<data_t>::getDataVector() const
const DataContainer<data_t>& LinearResidual<data_t>::getDataVector() const
{
if (!_dataVector)
throw std::runtime_error("LinearResidual::getDataVector: data vector not present");
......@@ -63,7 +63,6 @@ namespace elsa
return *_dataVector;
}
template <typename data_t>
LinearResidual<data_t>* LinearResidual<data_t>::cloneImpl() const
{
......
......@@ -59,7 +59,7 @@ namespace elsa
bool hasDataVector() const;
/// return the operator A (throws if the residual has none)
LinearOperator<data_t>& getOperator() const;
const LinearOperator<data_t>& getOperator() const;
/// return the data vector b (throws if the residual has none)
const DataContainer<data_t>& getDataVector() const;
......
......@@ -13,7 +13,7 @@ namespace elsa
*
* \tparam data_t data type for the domain of the residual of the functional, defaulting to real_t
*
* The Pseudohuber norm evaluates to \f$ \sum_{i=1}^n \delta \left( \sqrt{1 + (x_i / delta)^2} - 1 \right) \f$
* The Pseudohuber norm evaluates to \f$ \sum_{i=1}^n \delta \left( \sqrt{1 + (x_i / \delta)^2} - 1 \right) \f$
* for \f$ x=(x_i)_{i=1}^n \f$ and a slope parameter \f$ \delta \f$.
*
* Reference: https://doi.org/10.1109%2F83.551699
......
#include "Quadric.h"
#include "Identity.h"
#include <stdexcept>
......@@ -8,40 +9,81 @@ namespace elsa
Quadric<data_t>::Quadric(const LinearOperator<data_t>& A, const DataContainer<data_t>& b)
: Functional<data_t>(A.getDomainDescriptor()),
_linearResidual{A, b}
{}
template <typename data_t>
Quadric<data_t>::Quadric(const LinearOperator<data_t>& A)
: Functional<data_t>(A.getDomainDescriptor()),
_linearResidual{A}
{}
template <typename data_t>
Quadric<data_t>::Quadric(const DataContainer<data_t>& b)
: Functional<data_t>(b.getDataDescriptor()),
_linearResidual{b}
{}
template <typename data_t>
Quadric<data_t>::Quadric(const DataDescriptor& domainDescriptor)
: Functional<data_t>(domainDescriptor),
_linearResidual{domainDescriptor}
{}
template <typename data_t>
const LinearResidual<data_t>& Quadric<data_t>::getGradientExpression() const
{
// TODO: enable this spd check
// if (!A.isSpd())
// throw std::invalid_argument("Quadric: operator A needs to be spd");
return _linearResidual;
}
template <typename data_t>
data_t Quadric<data_t>::_evaluate(const DataContainer<data_t>& Rx)
{
auto temp = _linearResidual.getOperator().apply(Rx);
return static_cast<data_t>(0.5) * Rx.dot(temp) - Rx.dot(_linearResidual.getDataVector());
data_t xtAx;
if (_linearResidual.hasOperator()) {
auto temp = _linearResidual.getOperator().apply(Rx);
xtAx = Rx.dot(temp);
}
else {
xtAx = Rx.squaredL2Norm();
}
if (_linearResidual.hasDataVector()) {
return static_cast<data_t>(0.5) * xtAx - Rx.dot(_linearResidual.getDataVector());
}
else {
return static_cast<data_t>(0.5) * xtAx;
}
}
template <typename data_t>
void Quadric<data_t>::_getGradientInPlace(DataContainer<data_t>& Rx)
{
auto temp = _linearResidual.getOperator().apply(Rx);
temp -= _linearResidual.getDataVector();
Rx = temp;
Rx = _linearResidual.evaluate(Rx);
}
template <typename data_t>
LinearOperator<data_t> Quadric<data_t>::_getHessian(const DataContainer<data_t>& Rx)
{
return leaf(_linearResidual.getOperator());
if(_linearResidual.hasOperator())
return leaf(_linearResidual.getOperator());
else
return leaf(Identity<data_t>(*_domainDescriptor));
}
template <typename data_t>
Quadric<data_t>* Quadric<data_t>::cloneImpl() const
{
return new Quadric<data_t>(_linearResidual.getOperator(), _linearResidual.getDataVector());
if (_linearResidual.hasOperator() && _linearResidual.hasDataVector())
return new Quadric<data_t>(_linearResidual.getOperator(), _linearResidual.getDataVector());
else if (_linearResidual.hasOperator() && !_linearResidual.hasDataVector())
return new Quadric<data_t>(_linearResidual.getOperator());
else if (!_linearResidual.hasOperator() && _linearResidual.hasDataVector())
return new Quadric<data_t>(_linearResidual.getDataVector());
else
return new Quadric<data_t>(*_domainDescriptor);
}
template <typename data_t>
......
......@@ -11,10 +11,11 @@ namespace elsa
* \author Matthias Wieczorek - initial code
* \author Maximilian Hornung - modularization
* \author Tobias Lasser - modernization
* \author Nikola Dinev - add functionality
*
* \tparam data_t data type for the domain of the residual of the functional, defaulting to real_t
*
* The Quadric functional evaluates to \f$ \frac{1}{2} x^tAx + x^tb \f$ for a symmetric positive definite
* The Quadric functional evaluates to \f$ \frac{1}{2} x^tAx - x^tb \f$ for a symmetric positive definite
* operator A and a vector b.
*
* Please note: contrary to other functionals, Quadric does not allow wrapping an explicit residual.
......@@ -30,9 +31,33 @@ namespace elsa
*/
Quadric(const LinearOperator<data_t>& A, const DataContainer<data_t>& b);
/**
* \brief Constructor for the Quadric functional \f$ \frac{1}{2} x^tAx \f$ (trivial data vector)
*
* \param[in] A the operator (has to be symmetric positive definite)
*/
explicit Quadric(const LinearOperator<data_t>& A);
/**
* \brief Constructor for the Quadric functional \f$ \frac{1}{2} x^tx - x^tb \f$ (trivial operator)
*
* \param[in] b the data vector
*/
explicit Quadric(const DataContainer<data_t>& b);
/**
* \brief Constructor for the Quadric functional \f$ \frac{1}{2} x^tx \f$ (trivial operator and data vector)
*
* \param[in] domainDescriptor the descriptor of x
*/
explicit Quadric(const DataDescriptor& domainDescriptor);
/// default destructor
~Quadric() override = default;
/// returns the residual \f$ Ax - b \f$, which also corresponds to the gradient of the functional
const LinearResidual<data_t>& getGradientExpression() const;
protected:
/// the evaluation of the Quadric functional
data_t _evaluate(const DataContainer<data_t>& Rx) override;
......@@ -52,6 +77,9 @@ namespace elsa
private:
/// storing A,b in a linear residual
LinearResidual<data_t> _linearResidual;
/// lift from base class
using Functional<data_t>::_domainDescriptor;
};
} // namespace elsa
......@@ -8,13 +8,13 @@ namespace elsa
template <typename data_t>
WeightedL2NormPow2<data_t>::WeightedL2NormPow2(const Scaling<data_t>& weightingOp)
: Functional<data_t>(weightingOp.getDomainDescriptor()),
_weightingOp{weightingOp.clone()}
_weightingOp{static_cast<Scaling<data_t>*>(weightingOp.clone().release())}
{}
template <typename data_t>
WeightedL2NormPow2<data_t>::WeightedL2NormPow2(const Residual<data_t>& residual, const Scaling<data_t>& weightingOp)
: Functional<data_t>(residual),
_weightingOp{weightingOp.clone()}
_weightingOp{static_cast<Scaling<data_t>*>(weightingOp.clone().release())}
{
// sanity check
if (residual.getDomainDescriptor() != weightingOp.getDomainDescriptor() ||
......@@ -22,6 +22,12 @@ namespace elsa
throw std::invalid_argument("WeightedL2NormPow2: sizes of residual and weighting operator do not match");
}
template <typename data_t>
const Scaling<data_t>& WeightedL2NormPow2<data_t>::getWeightingOperator() const
{
return *_weightingOp;
}
template <typename data_t>
data_t WeightedL2NormPow2<data_t>::_evaluate(const DataContainer<data_t>& Rx)
{
......
......@@ -14,8 +14,8 @@ namespace elsa
*
* \tparam data_t data type for the domain of the functional, defaulting to real_t
*
* The weighted, squared l2 norm functional evaluates to \f$ 0.5 * \langle Wx, x \rangle \f$ using the
* standard scalar product, and where W is a diagonal scaling operator.
* The weighted, squared l2 norm functional evaluates to \f$ 0.5 * \| x \|_{W,2} = 0.5 * \langle x, Wx \rangle \f$
* using the standard scalar product, and where W is a diagonal scaling operator.
*/
template <typename data_t = real_t>
class WeightedL2NormPow2 : public Functional<data_t> {
......@@ -38,6 +38,9 @@ namespace elsa
/// default destructor
~WeightedL2NormPow2() override = default;
/// returns the weighting operator
const Scaling<data_t>& getWeightingOperator() const;
protected:
/// the evaluation of the weighted, squared l2 norm
data_t _evaluate(const DataContainer<data_t>& Rx) override;
......@@ -56,7 +59,7 @@ namespace elsa
private:
/// the weighting operator
std::unique_ptr<LinearOperator<data_t>> _weightingOp;
std::unique_ptr<Scaling<data_t>> _weightingOp;
};
} // namespace elsa
......@@ -21,11 +21,11 @@ public:
: LinearOperator<data_t>(domain, range) {}
protected:
void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) override {
void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override {
Ax = 1;
}
void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) override {
void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) const override {
Aty = 3;
}
......
......@@ -12,10 +12,137 @@
#include <catch2/catch.hpp>
#include "Quadric.h"
#include "Identity.h"
#include "Scaling.h"
using namespace elsa;
SCENARIO("Testing the Quadric functional") {
GIVEN("no operator and no data") {
IndexVector_t numCoeff(3); numCoeff << 13, 11, 7;
DataDescriptor dd(numCoeff);
WHEN("instantiating") {
Quadric func(dd);
THEN("the functional is as expected") {
REQUIRE(func.getDomainDescriptor() == dd);
auto *linRes = dynamic_cast<const LinearResidual<real_t> *>(&func.getResidual());
REQUIRE(linRes);
REQUIRE(linRes->hasDataVector() == false);
REQUIRE(linRes->hasOperator() == false);
const auto& gradExpr = func.getGradientExpression();
REQUIRE(gradExpr.hasDataVector() == false);
REQUIRE(gradExpr.hasOperator() == false);
}
THEN("a clone behaves as expected") {
auto qClone = func.clone();
REQUIRE(qClone.get() != &func);
REQUIRE(*qClone == func);
}
THEN("the evaluate, gradient and Hessian work as expected") {
RealVector_t dataVec(dd.getNumberOfCoefficients());
dataVec.setRandom();
DataContainer x(dd, dataVec);
real_t trueValue = static_cast<real_t>(0.5)*x.squaredL2Norm();
REQUIRE(func.evaluate(x) == Approx(trueValue));
REQUIRE(func.getGradient(x) == x);
REQUIRE(func.getHessian(x) == leaf(Identity(dd)));
}
}
}
GIVEN("an operator but no data") {
IndexVector_t numCoeff(3); numCoeff << 13, 11, 7;
DataDescriptor dd(numCoeff);
Scaling scalingOp(dd,static_cast<real_t>(3.0));
WHEN("instantiating") {
Quadric func(scalingOp);
THEN("the functional is as expected") {
REQUIRE(func.getDomainDescriptor() == dd);
auto *linRes = dynamic_cast<const LinearResidual<real_t> *>(&func.getResidual());
REQUIRE(linRes);
REQUIRE(linRes->hasDataVector() == false);
REQUIRE(linRes->hasOperator() == false);
const auto& gradExpr = func.getGradientExpression();
REQUIRE(gradExpr.hasDataVector() == false);
REQUIRE(gradExpr.getOperator() == scalingOp);
}
THEN("a clone behaves as expected") {
auto qClone = func.clone();
REQUIRE(qClone.get() != &func);
REQUIRE(*qClone == func);
}
THEN("the evaluate, gradient and Hessian work as expected") {
RealVector_t dataVec(dd.getNumberOfCoefficients());
dataVec.setRandom();
DataContainer x(dd, dataVec);
real_t trueValue = static_cast<real_t>(0.5)*scalingOp.getScaleFactor()*x.squaredL2Norm();
REQUIRE(func.evaluate(x) == Approx(trueValue));
REQUIRE(func.getGradient(x) == scalingOp.getScaleFactor()*x);
REQUIRE(func.getHessian(x) == leaf(scalingOp));
}
}
}
GIVEN("data but no operator") {
IndexVector_t numCoeff(3); numCoeff << 13, 11, 7;
DataDescriptor dd(numCoeff);
RealVector_t randomData(dd.getNumberOfCoefficients());
randomData.setRandom();
DataContainer dc(dd, randomData);
WHEN("instantiating") {
Quadric func(dc);
THEN("the functional is as expected") {
REQUIRE(func.getDomainDescriptor() == dd);
auto *linRes = dynamic_cast<const LinearResidual<real_t> *>(&func.getResidual());
REQUIRE(linRes);
REQUIRE(linRes->hasDataVector() == false);
REQUIRE(linRes->hasOperator() == false);
const auto& gradExpr = func.getGradientExpression();
REQUIRE(gradExpr.getDataVector() == dc);
REQUIRE(gradExpr.hasOperator() == false);
}
THEN("a clone behaves as expected") {
auto qClone = func.clone();
REQUIRE(qClone.get() != &func);
REQUIRE(*qClone == func);
}
THEN("the evaluate, gradient and Hessian work as expected") {
RealVector_t dataVec(dd.getNumberOfCoefficients());
dataVec.setRandom();
DataContainer x(dd, dataVec);
real_t trueValue = static_cast<real_t>(0.5)*x.squaredL2Norm() - x.dot(dc);
REQUIRE(func.evaluate(x) == Approx(trueValue));
REQUIRE(func.getGradient(x) == x - dc);
REQUIRE(func.getHessian(x) == leaf(Identity<real_t>(dd)));
}
}
}
GIVEN("an operator and data") {
IndexVector_t numCoeff(3); numCoeff << 13, 11, 7;
DataDescriptor dd(numCoeff);
......@@ -36,6 +163,10 @@ SCENARIO("Testing the Quadric functional") {
REQUIRE(linRes);
REQUIRE(linRes->hasDataVector() == false);
REQUIRE(linRes->hasOperator() == false);
const auto& gradExpr = func.getGradientExpression();
REQUIRE(gradExpr.getDataVector() == dc);
REQUIRE(gradExpr.getOperator() == idOp);
}
THEN("a clone behaves as expected") {
......
......@@ -31,6 +31,7 @@ SCENARIO("Testing the weighted, squared l2 norm functional") {
THEN("the functional is as expected") {
REQUIRE(func.getDomainDescriptor() == dd);
REQUIRE(func.getWeightingOperator() == scalingOp);
auto* linRes = dynamic_cast<const LinearResidual<real_t>*>(&func.getResidual());
REQUIRE(linRes);
......@@ -86,6 +87,7 @@ SCENARIO("Testing the weighted, squared l2 norm functional") {
THEN("the functional is as expected") {
REQUIRE(func.getDomainDescriptor() == dd);
REQUIRE(func.getWeightingOperator() == scalingOp);
auto* lRes = dynamic_cast<const LinearResidual<real_t>*>(&func.getResidual());