/usr/include/sopt/reweighted.h is in libsopt-dev 2.0.0-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 | #ifndef SOPT_REWEIGHTED_H
#define SOPT_REWEIGHTED_H
#include "sopt/linear_transform.h"
#include "sopt/types.h"
namespace sopt {
namespace algorithm {
template <class ALGORITHM> class Reweighted;
//! Factory function to create an l0-approximation by reweighting an l1 norm
template <class ALGORITHM>
Reweighted<ALGORITHM>
reweighted(ALGORITHM const &algo, typename Reweighted<ALGORITHM>::t_SetWeights const &set_weights,
typename Reweighted<ALGORITHM>::t_Reweightee const &reweightee);
//! \brief L0-approximation algorithm, through reweighting
//! \details This algorithm approximates \f$min_x ||Ψ^Tx||_0 + f(x)\f$ by solving the set of
//! problems \f$j\f$, \f$min_x ||W_jΨ^Tx||_1 + f(x)\f$ where the *diagonal* matrix \f$W_j\f$ is set
//! using the results from \f$j-1\f$: \f$ δ_j W_j^{-1} = δ_j + ||W_{j-1}Ψ^T||_1\f$. \f$δ_j\f$
//! prevents division by zero. It is a series which converges to zero. By default,
//! \f$δ_{j+1}=0.1δ_j\f$.
//!
//! The algorithm proceeds needs three forms of input:
//! - the inner algorithm, e.g. ImagingProximalADMM
//! - a function returning Ψ^Tx given x
//! - a function to modify the inner algorithm with new weights
template <class ALGORITHM> class Reweighted {
public:
//! Inner-loop algorithm
typedef ALGORITHM Algorithm;
//! Scalar type
typedef typename Algorithm::Scalar Scalar;
//! Real type
typedef typename real_type<Scalar>::type Real;
//! Weight vector type
typedef Vector<Real> WeightVector;
//! Type of then underlying vectors
typedef typename Algorithm::t_Vector XVector;
//! Type of the convergence function
typedef typename Algorithm::t_IsConverged t_IsConverged;
//! \brief Type of the function that is subject to reweighting
//! \details E.g. \f$Ψ^Tx\f$. Note that l1-norm is not applied here.
typedef std::function<XVector(Algorithm const &, XVector const &)> t_Reweightee;
//! Type of the function to set weights
typedef std::function<void(Algorithm &, WeightVector const &)> t_SetWeights;
//! Function to update delta at each turn
typedef std::function<Real(Real)> t_DeltaUpdate;
//! output from running reweighting scheme
struct ReweightedResult {
//! Number of iterations (outer loop)
t_uint niters;
//! Wether convergence was achieved
bool good;
//! Weights at last iteration
WeightVector weights;
//! Result from last inner loop
typename Algorithm::DiagnosticAndResult algo;
//! Default construction
ReweightedResult() : niters(0), good(false), weights(WeightVector::Ones(1)), algo() {}
};
Reweighted(Algorithm const &algo, t_SetWeights const &setweights, t_Reweightee const &reweightee)
: algo_(algo), setweights_(setweights), reweightee_(reweightee),
itermax_(std::numeric_limits<t_uint>::max()), min_delta_(0e0), is_converged_(),
update_delta_([](Real delta) { return 1e-1 * delta; }) {}
//! Underlying "inner-loop" algorithm
Algorithm &algorithm() { return algo_; }
//! Underlying "inner-loop" algorithm
Algorithm const &algorithm() const { return algo_; }
//! Sets the underlying "inner-loop" algorithm
Reweighted<Algorithm> &algorithm(Algorithm const &algo) {
algo_ = algo;
return *this;
}
//! Sets the underlying "inner-loop" algorithm
Reweighted<Algorithm> &algorithm(Algorithm &&algo) {
algo_ = std::move(algo);
return *this;
}
//! Function to reset the weights in the algorithm
t_SetWeights const &set_weights() const { return setweights_; }
//! Function to reset the weights in the algorithm
Reweighted<Algorithm> &set_weights(t_SetWeights const &setweights) const {
setweights_ = setweights;
return *this;
}
//! Sets the weights on the underlying algorithm
void set_weights(Algorithm &algo, WeightVector const &weights) const {
return set_weights()(algo, weights);
}
//! Function that needs to be reweighted
//! \details E.g. \f$Ψ^Tx\f$. Note that l1-norm is not applied here.
Reweighted<Algorithm> &reweightee(t_Reweightee const &rw) {
reweightee_ = rw;
return *this;
}
//! Function that needs to be reweighted
t_Reweightee const &reweightee() const { return reweightee_; }
//! Forwards to the reweightee function
XVector reweightee(XVector const &x) const { return reweightee()(algorithm(), x); }
//! Maximum number of reweighted iterations
t_uint itermax() const { return itermax_; }
Reweighted &itermax(t_uint i) {
itermax_ = i;
return *this;
}
//! Lower limit for delta
Real min_delta() const { return min_delta_; }
Reweighted &min_delta(Real min_delta) {
min_delta_ = min_delta;
return *this;
}
//! Checks convergence of the reweighting scheme
t_IsConverged const &is_converged() const { return is_converged_; }
Reweighted &is_converged(t_IsConverged const &convergence) {
is_converged_ = convergence;
return *this;
}
bool is_converged(XVector const &x) const { return is_converged() ? is_converged()(x) : false; }
//! \brief Performs reweighting
//! \details This overload will compute an initial result without initial weights set to one.
template <class INPUT>
typename std::enable_if<not(std::is_same<INPUT, typename Algorithm::DiagnosticAndResult>::value
or std::is_same<INPUT, ReweightedResult>::value),
ReweightedResult>::type
operator()(INPUT const &input) const;
//! \brief Performs reweighting
//! \details This overload will compute an initial result without initial weights set to one.
ReweightedResult operator()() const;
//! Reweighted algorithm, from prior call to inner-algorithm
ReweightedResult operator()(typename Algorithm::DiagnosticAndResult const &warm) const;
//! Reweighted algorithm, from prior call to reweighting algorithm
ReweightedResult operator()(ReweightedResult const &warm) const;
//! Updates delta
Real update_delta(Real delta) const { return update_delta()(delta); }
//! Updates delta
t_DeltaUpdate const &update_delta() const { return update_delta_; }
//! Updates delta
Reweighted<Algorithm> update_delta(t_DeltaUpdate const &ud) const { return update_delta_ = ud; }
protected:
//! Inner loop algorithm
Algorithm algo_;
//! Function to set weights
t_SetWeights setweights_;
//! \brief Function that is subject to reweighting
//! \details E.g. \f$Ψ^Tx\f$. Note that l1-norm is not applied here.
t_Reweightee reweightee_;
//! Maximum number of reweighted iterations
t_uint itermax_;
//! \brief Lower limit for delta
Real min_delta_;
//! Checks convergence
t_IsConverged is_converged_;
//! Updates delta at each turn
t_DeltaUpdate update_delta_;
};
template <class ALGORITHM>
template <class INPUT>
typename std::
enable_if<not(std::is_same<INPUT, typename ALGORITHM::DiagnosticAndResult>::value
or std::is_same<INPUT, typename Reweighted<ALGORITHM>::ReweightedResult>::value),
typename Reweighted<ALGORITHM>::ReweightedResult>::type
Reweighted<ALGORITHM>::operator()(INPUT const &input) const {
Algorithm algo = algorithm();
set_weights(algo, WeightVector::Ones(1));
return operator()(algo(input));
}
template <class ALGORITHM>
typename Reweighted<ALGORITHM>::ReweightedResult Reweighted<ALGORITHM>::operator()() const {
Algorithm algo = algorithm();
set_weights(algo, WeightVector::Ones(1));
return operator()(algo());
}
template <class ALGORITHM>
typename Reweighted<ALGORITHM>::ReweightedResult Reweighted<ALGORITHM>::
operator()(typename Algorithm::DiagnosticAndResult const &warm) const {
ReweightedResult result;
result.algo = warm;
result.weights = WeightVector::Ones(1);
return operator()(result);
}
template <class ALGORITHM>
typename Reweighted<ALGORITHM>::ReweightedResult Reweighted<ALGORITHM>::
operator()(ReweightedResult const &warm) const {
SOPT_HIGH_LOG("Starting reweighted scheme");
// Copies inner algorithm, so that operator() can be constant
Algorithm algo(algorithm());
ReweightedResult result(warm);
auto delta = std::max(standard_deviation(reweightee(warm.algo.x)), min_delta());
SOPT_LOW_LOG("- Initial delta: {}", delta);
for(result.niters = 0; result.niters < itermax(); ++result.niters) {
SOPT_LOW_LOG("Reweigting iteration {}/{} ", result.niters, itermax());
SOPT_LOW_LOG(" - delta: {}", delta);
result.weights = delta / (delta + reweightee(result.algo.x).array().abs());
set_weights(algo, result.weights);
result.algo = algo(result.algo);
if(is_converged(result.algo.x)) {
SOPT_MEDIUM_LOG("Reweighting scheme did converge in {} iterations", result.niters);
result.good = true;
break;
}
delta = std::max(min_delta(), update_delta(delta));
}
// result is always good if no convergence function is defined
if(not is_converged())
result.good = true;
else if(not result.good)
SOPT_ERROR("Reweighting scheme did *not* converge in {} iterations", itermax());
return result;
}
//! Factory function to create an l0-approximation by reweighting an l1 norm
template <class ALGORITHM>
Reweighted<ALGORITHM>
reweighted(ALGORITHM const &algo, typename Reweighted<ALGORITHM>::t_SetWeights const &set_weights,
typename Reweighted<ALGORITHM>::t_Reweightee const &reweightee) {
return {algo, set_weights, reweightee};
}
template <class SCALAR> class ImagingProximalADMM;
template <class ALGORITHM> class PositiveQuadrant;
template <class T>
Eigen::CwiseUnaryOp<const details::ProjectPositiveQuadrant<typename T::Scalar>, const T>
positive_quadrant(Eigen::DenseBase<T> const &input);
template <class SCALAR>
Reweighted<PositiveQuadrant<ImagingProximalADMM<SCALAR>>>
reweighted(ImagingProximalADMM<SCALAR> const &algo) {
auto const posq = positive_quadrant(algo);
typedef typename std::remove_const<decltype(posq)>::type Algorithm;
typedef Reweighted<Algorithm> RW;
auto const reweightee
= [](Algorithm const &posq, typename RW::XVector const &x) -> typename RW::XVector {
return posq.algorithm().Psi().adjoint() * x;
};
auto const set_weights = [](Algorithm &posq, typename RW::WeightVector const &weights) -> void {
posq.algorithm().l1_proximal_weights(weights);
};
return {posq, set_weights, reweightee};
}
} // namespace algorithm
} // namespace sopt
#endif
|