atoms.hpp

Go to the documentation of this file.

#ifndef MLPACK_CORE_OPTIMIZERS_FW_ATOMS_HPP
#define MLPACK_CORE_OPTIMIZERS_FW_ATOMS_HPP

#include <mlpack/prereqs.hpp>
#include <mlpack/core/optimizers/proximal/proximal.hpp>
#include "func_sq.hpp"

namespace mlpack {
namespace optimization {

class Atoms
{
 public:
  Atoms(){ /* Nothing to do. */ }

  void AddAtom(const arma::vec& v, FuncSq& function, const double c = 0)
  {
    if (currentAtoms.is_empty())
    {
      CurrentAtoms() = v;
      CurrentCoeffs().set_size(1);
      CurrentCoeffs().fill(c);
      atomSqTerm.set_size(1);
      atomSqTerm(0) = std::pow(norm(function.MatrixA() * v, 2), 2);
    }
    else
    {
      currentAtoms.insert_cols(0, v);
      arma::vec cVec(1);
      cVec(0) = c;
      currentCoeffs.insert_rows(0, cVec);
      double tmp = std::pow(norm(function.MatrixA() * v, 2), 2);
      arma::vec tmpVec(1);
      tmpVec(0) = tmp;
      atomSqTerm.insert_rows(0, tmpVec);
    }
  }


  void RecoverVector(arma::mat& x)
  {
    x = currentAtoms * currentCoeffs;
  }

  void PruneSupport(const double F, FuncSq& function)
  {
    arma::vec sqTerm = 0.5 * atomSqTerm % square(currentCoeffs);

    while (currentAtoms.n_cols > 1)
    {
      // Solve for current gradient.
      arma::mat x;
      RecoverVector(x);
      arma::mat gradient(size(x));
      function.Gradient(x, gradient);

      // Find possible atom to be deleted.
      arma::vec gap = sqTerm -
          currentCoeffs % trans(gradient.t() * currentAtoms);
      arma::uword ind;
      gap.min(ind);

      // Try deleting the atom.
      arma::mat newAtoms(currentAtoms.n_rows, currentAtoms.n_cols - 1);
      if (ind > 0)
        newAtoms.cols(0, ind - 1) = currentAtoms.cols(0, ind - 1);
      if (ind < (currentAtoms.n_cols - 1))
      {
        newAtoms.cols(ind, newAtoms.n_cols - 1) =
            currentAtoms.cols(ind + 1, currentAtoms.n_cols - 1);
      }

      // Reoptimize the coefficients, we brute-forcely reoptimize in the span,
      // which would be used in UpdateSpan class. Alternatively, if you want to
      // add an atom norm constraint, you could use projected gradient method,
      // see the implementaton of ProjectedGradientEnhancement().
      arma::vec newCoeffs =
          solve(function.MatrixA() * newAtoms, function.Vectorb());

      // Evaluate the function again.
      double Fnew = function.Evaluate(newAtoms * newCoeffs);

      if (Fnew > F)
        // Should not delete the atom.
        break;
      else
      {
        // Delete the atom from current atoms.
        currentAtoms = newAtoms;
        currentCoeffs = newCoeffs;
        atomSqTerm.shed_row(ind);
        sqTerm.shed_row(ind);
      } // else
    } // while
  }


  void ProjectedGradientEnhancement(FuncSq& function,
                                    double tau,
                                    double stepSize,
                                    size_t maxIteration = 100,
                                    double tolerance = 1e-3)
  {
    arma::mat x;
    RecoverVector(x);
    double value = function.Evaluate(x);

    for (size_t iter = 1; iter<maxIteration; iter++)
    {
      // Update currentCoeffs with gradient descent method.
      arma::mat g;
      function.Gradient(x, g);
      g = currentAtoms.t() * g;
      currentCoeffs = currentCoeffs - stepSize * g;

      // Projection of currentCoeffs to satisfy the atom norm constraint.
      Proximal::ProjectToL1Ball(currentCoeffs, tau);

      RecoverVector(x);
      double valueNew = function.Evaluate(x);

      if ((value - valueNew) < tolerance)
        break;

      value = valueNew;
    }
  }


  const arma::vec& CurrentCoeffs() const { return currentCoeffs; }
  arma::vec& CurrentCoeffs() { return currentCoeffs; }

  const arma::mat& CurrentAtoms() const { return currentAtoms; }
  arma::mat& CurrentAtoms() { return currentAtoms; }

 private:
  arma::vec currentCoeffs;

  arma::mat currentAtoms;

  arma::vec atomSqTerm;
}; // class Atoms
}  // namespace optimization
}  // namespace mlpack

#endif