lostella / libforbes Goto Github PK

Implement the indicator of the projection onto the 1-norm using the code of J. Duchi:

function w = ProjectOntoL1Ball(v, b)
if (norm(v, 1) < b)
  w = v;
  return;
end
u = sort(abs(v),'descend');
sv = cumsum(u);
rho = find(u > (sv - b) ./ (1:length(u))', 1, 'last');
theta = max(0, (sv(rho) - b) / rho);
w = sign(v) .* max(abs(v) - theta, 0);

LeastSquares.cpp is missing. Also, it is listed in the wrong section of the Makefile I believe.

As a side note: I kind of guess what this function is like, for the moment I'm using QuadraticLoss for implementing least squares problems (in tests).

We need to abolish the method virtual int call(Matrix& x, double& f, Matrix& grad, Matrix& hessian); and replace it by virtual int call(Matrix& x, double& f, Matrix& grad, LinearOperator& hessian);, where the Hessian at x is returned as an instance of LinearOperator

TODOs:

1. Remove virtual int call(Matrix& x, double& f, Matrix& grad, Matrix& hessian); in Function.h
2. Replace it with virtual int call(Matrix& x, double& f, Matrix& grad, LinearOperator& hessian);
3. Modify the corresponding method in Quadratic

Implement the max{x,0} operator in Matrix like

   void plusop();
   void plusop(Matrix * mat);

where the former will update the current matrix x with [x]+, and the latter will update a given matrix mat with the result [x]+.

Method int FBSplitting::setMaxIt(int) does not return anything (see FBSplitting.cpp:26).

We need a constructor in IterativeSolver which needs to properly initialise m_it and m_maxit. Derived classes, such as FBSplitting should extend this constructor or properly override it.

Also, document the default tolerance.

Why do we store m_res1_rows and m_res1_cols as private fields in FBCache? Since we store m_L1 we can retrieve them from there any time we like. The same holds for m_x (we store m_x_rows and m_x_cols).

Further, we don't even need to store m_L1, m_f1, etc. This is why you introduced FBProblem as a field in FBCache. We have the constructor FBCache(FBProblem & p, Matrix & x, double gamma); so we can retrieve a reference to an FBProblem and store it internally and this will be it.

I believe we can/should refactor FBCache and simplify it.

Meticulously check returned status codes. For instance, in IterativeSolver we now have

int IterativeSolver::run() {
    int status = ForBESUtils::STATUS_OK;
    while (m_it < m_maxit && !stop() && !ForBESUtils::is_status_error(status)) {
        status = iterate();
        m_it++;
    }
    return status;
}

We should take care of status codes in FBCache (see this cppcheck report for FBCache). We should add assertions in tests to check whether the status codes returned from function invocations and solvers are STATUS_OK.

Implement and test the indicator of an affine subspace defined as Ax=b. For that, we need to solve a least-squares problem which may require to implement a QR decomposition class or a least squares solver.

Implement virtual int call(Matrix& x, double& f, Matrix& grad, LinearOperator& hessian); in Quadratic and LogLogisticLoss. For the former, this should be a minor amendment.

Occasionally, the following tests fail (no idea why!). By the way, this is the only test in the test suite that fails. Can you confirm this? This started happening after beffc81.

Test name: TestFBSplitting::testBoxQP_small
assertion failed
- Expression: solver->getIt() < max_iter

TestFBSplitting.cpp:127:Assertion
Test name: TestFBSplitting::testLasso_small
assertion failed
- Expression: solver->getIt() < max_iter

TestFBSplitting.cpp:173:Assertion
Test name: TestFBSplitting::testSparseLogReg_small
assertion failed
- Expression: solver->getIt() < max_iter

Failures !!!
Run: 3   Failure total: 3   Failures: 3   Errors: 0
make: *** [test] Error 1

In FBCache private members m_FBEx and m_sqnormFPRx are not initialised in the constructor.

Method int IterativeSolver::run() does not return anything in IterativeSolver.cpp:9.

Method int FBSplitting::setTol(double) does not return anything (see FBSplitting.cpp:30)

The tolerance parameter, m_tol, in FBSplitting is not initialised in the constructor. A default value should be assigned to it. I believe we should have an additional constructor to specify the tolerance directly. What is more, it seems to me that m_maxit and m_tol could be moved to IterativeSolver.

When the function QuadraticOverAffine(Q, q, A, b) is defined, then its conjugate is quadratic as long as d'Qd > 0 for all d in the subspace {x : Ax = b}. Now when solving

    minimize f(x) + g(z) subject to Ax + Bz = b

we will assume that f is strongly convex (this will be our contract with the user, so to speak) and therefore if f = QuadraticOverAffine(Q, q, A, b) then necessarily f* will be quadratic. This should be easily detectable: just like we mark whether or not a function f defines the gradient of the conjugate and then set that ConjugateFunction(f) defines the gradient when constructing it, it would be desirable to mark whether f is "conjugate-quadratic" and set that ConjugateFunction(f) is quadratic.

Implement and test IndProbSimplex (both call and callProx).

Script install.sh is omits the installation of libcppunit, so users who do a fresh installation will not be able to run the tests.

The following snippet runs an FBS algorithm for a QP problem with inequality and equality constraints. The primal solution is correct (I cross-checked with Gurobi), but the algorithm reaches the maximum number of iterations always (even for large tolerance values).

const size_t n = 4;
    double data_Q[] = {
        3, 1, 1, 0,
        1, 3, 0, -1,
        1, 0, 3, 0,
        0, -1, 0, 5
    };
    double data_q[] = {1, 2, 3, 4};
    Matrix Q(n, n, data_Q);
    Matrix q(n, 1, data_q);


    const size_t p = 2;
    double data_A[] = {
        1, 0,
        1, 0,
        0, 1,
        0, 2
    };
    double data_b[] = {
        1,
        0
    };

    Matrix A(p, n, data_A);
    Matrix b(p, 1, data_b);

    std::cout << Q;
    std::cout << q;
    std::cout << A;
    std::cout << b;

    // Define f and f*
    QuadOverAffine f(Q, q, A, b);
    ConjugateFunction f_conj(f);

    // Define g and g*
    double lb = -2;
    double ub = +2;
    IndBox g(lb, ub);
    ConjugateFunction g_conj(g);
    double data_H[] = {1, 1, 1, 1};
    Matrix H(1, n, data_H);
    H.transpose();
    H *= -1.0;
    MatrixOperator Op_Htr_minus(H);

    FBProblem prob(f_conj, Op_Htr_minus, g_conj);
    Matrix y0(1, 1);
    FBStoppingRelative stop(1e-3);
    FBSplitting fbs(prob, y0, 0.01, stop, 45000);
    fbs.run();

    Matrix ystar = fbs.getSolution();
    Matrix H_ystar = H*ystar;
    double f_conj_star;
    Matrix xstar(n, 1);
    f_conj.call(H_ystar, f_conj_star, xstar);

    std::cout << fbs.getIt();
    std::cout << xstar;

Implement CongruentSeparableSum, a fast variant of SeparableSum when the partitioning of the argument of f(x) is done in a congruent fashion.

In FBCache the following piece of code has the following problem:

if (m_prob.f1() != NULL) {
        Matrix gradf1_diff = Matrix(m_x_rows, m_x_cols);
        if (m_prob.L1() != NULL) {
            Matrix gradf1_L1diff = Matrix(m_prob.L1()->dimensionOut());
            gradf1_diff = m_prob.L1()->callAdjoint(gradf1_L1diff);
        } else {

        }
        Matrix::add(*m_gradFBEx, -1.0, gradf1_diff, 1.0 / gamma);
    }

The line

Matrix gradf1_diff = Matrix(m_x_rows, m_x_cols);

allocates memory for a dense matrix. Then, line

gradf1_diff = m_prob.L1()->callAdjoint(gradf1_L1diff);

computes the value of the adjoint of L1 at gradf1_L1diff and copies the value to gradf1_diff!
This happens with the invocation of Matrix::operator=. Notice that would be different to do

Matrix gradf1_diff = m_prob.L1()->callAdjoint(gradf1_L1diff);

Maybe we could move int setMaxIt(int maxit); from FBSplitting into IterativeSolver. All solvers define a maximum number of iterations, right?

The following Functions need to be implemented:

1. IndPos
2. IndNeg
3. IndAffineSubspace
4. IndProbSimplex
5. DistanceToBall2
6. CongruentSeparableSum

Implement new methods in Matrix to return the 2-norm/Frobenius norm and the maximum-absolute-value norm of vectors/matrices.