Major code reorganization

@FBCache/Check_Gamma.m → +forbes/+fbe/@FBCache/Check_Gamma.m

@@ -18,7 +18,7 @@
 function [flag, cache_pg] = Check_Gamma(cache, bet)
 
 z = cache.Get_ProxGradStep();
-cache_pg = FBCache(cache.prob, z, cache.gam, cache.ops);
+cache_pg = forbes.fbe.FBCache(cache.prob, z, cache.gam, cache.ops);
 fz = cache_pg.Get_f();
 
 % We should do the following, but maybe accessing the fields directly has

@FBCache/Get_CacheLine.m → +forbes/+fbe/@FBCache/Get_CacheLine.m

@@ -16,7 +16,7 @@
 gam = cache.gam;
 
 if nargin < 4
-    cachet = FBCache(cache.prob, cache.x + tau*cache.dir1, cache.gam, cache.ops);
+    cachet = forbes.fbe.FBCache(cache.prob, cache.x + tau*cache.dir1, cache.gam, cache.ops);
 end
 
 if nargin < 4 || ~cachet.flagGradStep
@@ -33,9 +33,9 @@
     if prob.istheref2
         cachet.res2x = cache.res2x + tau*cache.C2dir1;
         if prob.useHessian
-            [f2xt, gradf2res2xt, cachet.Hessf2res2x] = prob.callf2(cachet.res2x);
+            [gradf2res2xt, f2xt, cachet.Hessf2res2x] = prob.f2.gradient(cachet.res2x);
         else
-            [f2xt, gradf2res2xt] = prob.callf2(cachet.res2x);
+            [gradf2res2xt, f2xt] = prob.f2.gradient(cachet.res2x);
             cachet.gradf2res2x = gradf2res2xt;
         end
         if cache.flagOps
@@ -68,10 +68,10 @@
 if nargin < 4 || ~cachet.flagProxGradStep
     if prob.isthereD
         mugam = prob.mu*gam;
-        [z, cachet.gz] = prob.callg(prob.D*cachet.y, mugam);
+        [z, cachet.gz] = prob.g.prox(prob.D*cachet.y, mugam);
         cachet.z = cachet.y + prob.D'*(z - prob.D*cachet.y)/prob.mu;
     else
-        [cachet.z, cachet.gz] = prob.callg(cachet.y, gam);
+        [cachet.z, cachet.gz] = prob.g.prox(cachet.y, gam);
     end
     if cache.flagOps
         cache.ops.addproxg();

@FBCache/Get_CacheSegment.m → +forbes/+fbe/@FBCache/Get_CacheSegment.m

@@ -7,7 +7,7 @@
 prob = cache.prob;
 gam = cache.gam;
 
-cachet = FBCache(cache.prob, cache.x + tau*cache.dir1 + (1-tau)*cache.dir2, cache.gam, cache.ops);
+cachet = forbes.fbe.FBCache(cache.prob, cache.x + tau*cache.dir1 + (1-tau)*cache.dir2, cache.gam, cache.ops);
 
 fxt = 0;
 gradfxt = 0;
@@ -22,9 +22,9 @@
 if prob.istheref2
     cachet.res2x = cache.res2x + tau*cache.C2dir1 + (1-tau)*cache.C2dir2;
     if prob.useHessian
-        [f2xt, gradf2res2xt, cachet.Hessf2res2x] = prob.callf2(cachet.res2x);
+        [gradf2res2xt, f2xt, cachet.Hessf2res2x] = prob.f2.gradient(cachet.res2x);
     else
-        [f2xt, gradf2res2xt] = prob.callf2(cachet.res2x);
+        [gradf2res2xt, f2xt] = prob.f2.gradient(cachet.res2x);
         cachet.gradf2res2x = gradf2res2xt;
     end
     if cache.flagOps
@@ -54,10 +54,10 @@
 
 if prob.isthereD
     mugam = prob.mu*gam;
-    [z, cachet.gz] = prob.callg(prob.D*cachet.y, mugam);
+    [z, cachet.gz] = prob.g.prox(prob.D*cachet.y, mugam);
     cachet.z = cachet.y + prob.D'*(z - prob.D*cachet.y)/prob.mu;
 else
-    [cachet.z, cachet.gz] = prob.callg(cachet.y, gam);
+    [cachet.z, cachet.gz] = prob.g.prox(cachet.y, gam);
 end
 if cache.flagOps
     cache.ops.addproxg();

@FBCache/Get_GradFBE.m → +forbes/+fbe/@FBCache/Get_GradFBE.m

@@ -45,7 +45,7 @@
     else
         % imaginary trick
         res2xepsFPR = cache.res2x + 1e-100i*C2FPR;
-        [~, gradf2res2xepsd] = prob.callf2(res2xepsFPR);
+        [gradf2res2xepsd, ~] = prob.f2.gradient(res2xepsFPR);
         if cache.flagOps, cache.ops.addgradf2(); end
         HC2FPR = imag(gradf2res2xepsd)/1e-100;
         % forward differences

@FBCache/Get_Gradf.m → +forbes/+fbe/@FBCache/Get_Gradf.m

@@ -28,18 +28,18 @@
 if prob.istheref2
     if prob.isthereC2
         if prob.useHessian
-            [~, gradf2res2x, cache.Hessf2res2x] = prob.callf2(cache.res2x);
+            [gradf2res2x, ~, cache.Hessf2res2x] = prob.f2.gradient(cache.res2x);
         else
-            [~, gradf2res2x] = prob.callf2(cache.res2x);
+            [gradf2res2x, ~] = prob.f2.gradient(cache.res2x);
             cache.gradf2res2x = gradf2res2x;
         end
         cache.gradf2x = prob.C2'*gradf2res2x;
         if cache.flagOps, cache.ops.addC2(); end
     else
         if prob.useHessian
-            [~, gradf2res2x, cache.Hessf2res2x] = prob.callf2(cache.res2x);
+            [gradf2res2x, ~, cache.Hessf2res2x] = prob.f2.gradient(cache.res2x);
         else
-            [~, gradf2res2x] = prob.callf2(cache.res2x);
+            [gradf2res2x, ~] = prob.f2.gradient(cache.res2x);
             cache.gradf2res2x = gradf2res2x;
         end
         cache.gradf2x = gradf2res2x;

@FBCache/Get_ProxGradStep.m → +forbes/+fbe/@FBCache/Get_ProxGradStep.m

@@ -20,10 +20,10 @@
 
 if prob.isthereD
     mugam = prob.mu*gam;
-    [z, cache.gz] = prob.callg(prob.D*cache.y, mugam);
+    [z, cache.gz] = prob.g.prox(prob.D*cache.y, mugam);
     cache.z = cache.y + prob.D'*(z - prob.D*cache.y)/prob.mu;
 else
-    [cache.z, cache.gz] = prob.callg(cache.y, gam);
+    [cache.z, cache.gz] = prob.g.prox(cache.y, gam);
 end
 if cache.flagOps
     cache.ops.addproxg();

@FBCache/Get_f.m → +forbes/+fbe/@FBCache/Get_f.m

@@ -17,10 +17,10 @@
         C1x = prob.C1*cache.x;
         if cache.flagOps, cache.ops.addC1(); end
         cache.res1x = C1x + prob.d1;
-        [cache.f1x, cache.gradf1res1x] = prob.callf1(cache.res1x);
+        [cache.gradf1res1x, cache.f1x] = prob.f1.gradient(cache.res1x);
     else
         cache.res1x = cache.x + prob.d1;
-        [cache.f1x, cache.gradf1res1x] = prob.callf1(cache.res1x);
+        [cache.gradf1res1x, cache.f1x] = prob.f1.gradient(cache.res1x);
         cache.gradf1x = cache.gradf1res1x;
     end
     if cache.flagOps
@@ -36,10 +36,10 @@
         C2x = prob.C2*cache.x;
         if cache.flagOps, cache.ops.addC2(); end
         cache.res2x = C2x + prob.d2;
-        f2x = prob.callf2(cache.res2x);
+        [~, f2x] = prob.f2.gradient(cache.res2x);
     else
         cache.res2x = cache.x + prob.d2;
-        f2x = prob.callf2(cache.res2x);
+        [~, f2x] = prob.f2.gradient(cache.res2x);
     end
     if cache.flagOps, cache.ops.addf2(); end
     cache.f2x = f2x;

+forbes/+functions/@Conjugate/Conjugate.m

@@ -0,0 +1,27 @@
+% CONJUGATE Fenchel conjugate function (of some other given function)
+
+classdef Conjugate < forbes.functions.Proximable
+    properties
+        f
+    end
+    methods
+        function obj = Conjugate(f)
+            obj.f = f;
+        end
+        function p = is_convex(obj)
+            p = true;
+        end
+        function p = is_strongly_convex(obj)
+            p = obj.f.is_smooth();
+        end
+        function p = is_smooth(obj)
+            p = obj.f.is_strongly_convex();
+        end
+        function p = is_quadratic(obj)
+            p = obj.f.is_strongly_convex() && obj.f.is_generalized_quadratic();
+        end
+        function p = is_generalized_quadratic(obj)
+            p = obj.f.is_quadratic();
+        end
+    end
+end

+forbes/+functions/@Conjugate/compute_gradient.m

@@ -0,0 +1,3 @@
+function [g, v] = compute_gradient(obj, x)
+    [g, v] = obj.f.compute_gradient_conjugate(x);
+end

+forbes/+functions/@Conjugate/compute_prox.m

@@ -0,0 +1,8 @@
+function [p, v] = compute_prox(obj, x, gam)
+    % use Moreau identity:
+    %   prox(x; gam, f*) = x - gam*prox(x/gam; 1/gam, f)
+    %
+    [p1, v1] = obj.f.prox(x/gam, 1/gam);
+    p = x - gam*p1;
+    v = x(:)'*p1(:) - gam*p1(:)'*p1(:) - v1;
+end

+forbes/+functions/@DistBoxL1/DistBoxL1.m

@@ -0,0 +1,19 @@
+% DISTBOXL1 L1 distance from a box
+
+classdef DistBoxL1 < forbes.functions.Proximable
+    properties
+        lb, ub % box bounds
+        w
+    end
+    methods
+        function obj = DistBoxL1(lb, ub, w)
+            if nargin < 3, w = 1.0; end
+            obj.lb = lb;
+            obj.ub = ub;
+            obj.w = w;
+        end
+        function p = is_convex(obj)
+            p = true;
+        end
+    end
+end

+forbes/+functions/@DistBoxL1/compute_prox.m

@@ -0,0 +1,13 @@
+function [p, v] = compute_prox(obj, x, gam)
+    mu = gam * obj.w;
+    p = max(x-obj.ub-mu, 0) - max(obj.lb-x-mu, 0) + min(max(x, obj.lb), obj.ub);
+    if nargout > 1
+        proj = max(obj.lb, min(obj.ub, p));
+        if isscalar(obj.w)
+            v = sum(obj.w*abs(p-proj));
+        else
+            finw = ~isinf(obj.w);
+            v = sum(obj.w(finw).*abs(p(finw)-proj(finw)));
+        end
+    end
+end

+forbes/+functions/@DistL2/DistL2.m

@@ -0,0 +1,18 @@
+% DISTL2 L2 (Euclidean) distance from a convex set
+
+classdef DistL2 < forbes.functions.Proximable
+    properties
+        ind % indicator function of the (convex) set
+        lam
+    end
+    methods
+        function obj = DistL2(ind, lam)
+            if nargin < 2, lam = 1.0; end
+            obj.ind = ind;
+            obj.lam = lam;
+        end
+        function p = is_convex(obj)
+            p = true;
+        end
+    end
+end

+forbes/+functions/@DistL2/compute_prox.m

@@ -0,0 +1,13 @@
+function [y, v] = compute_prox(obj, x, gam)
+    p = obj.ind.compute_prox(x, 1.0);
+    d = norm(x-p, 'fro');
+    gamlam = (gam*obj.lam);
+    if gamlam < d
+        gamlamd = gamlam/d;
+        y = (1-gamlamd)*x + gamlamd*p;
+        v = obj.lam*(d-gamlam);
+    else
+        y = p;
+        v = 0.0;
+    end
+end

+forbes/+functions/@EpiCompose/EpiCompose.m

@@ -0,0 +1,33 @@
+% EPICOMPOSE Epi-composition
+%
+%   EPICOMPOSE(f, A) returns the function g(y) = inf {f(x) : Ax = y}.
+
+classdef EpiCompose < forbes.functions.Proximable
+    properties
+        f, A
+        mu
+        flag % f is 1: Quadratic, 2: QuadraticOverAffine, 3: Anything else
+        Q, q, C, d % to store data of (generalized) quadratic functions
+        gam_prox, L_prox % to store Cholesky factor
+        x
+    end
+    methods
+        function obj = EpiCompose(f, A)
+            if isa(f, 'Quadratic')
+                obj.flag = 1;
+                obj.Q = f.Q;
+                obj.q = f.q;
+                obj.gam_prox = 0;
+                if isscalar(A), A = A*speye(size(f.Q, 2)); end
+            else
+                obj.flag = 3;
+                obj.mu = Proximable.get_gram_diagonal(A');
+                if obj.mu == 0
+                    error('A must be s.t. A''*A = mu*Id');
+                end
+            end
+            obj.f = f;
+            obj.A = A;
+        end
+    end
+end

+forbes/+functions/@EpiCompose/compute_prox.m

@@ -0,0 +1,19 @@
+function [p, v] = compute_prox(obj, s, gam)
+    switch obj.flag
+        case 1 % f = Quadratic
+            if gam ~= obj.gam_prox % more robust test?
+                % factor matrix when gam changes (or at first call)
+                obj.gam_prox = gam;
+                obj.L_prox = chol(obj.Q + (1/gam)*obj.A'*obj.A); % do differently for sparse?
+            end
+            obj.x = obj.L_prox\(obj.L_prox'\((obj.A'*s)/gam - obj.q));
+            v = 0.5*(obj.x'*(obj.Q*obj.x)) + obj.q'*obj.x; % can we save something here?
+        case 2 % f = QuadraticOverAffine
+            error('not implemented');
+        case 3 % f = Any proximable function and A'*A = mu*Id, mu > 0
+            [obj.x, v] = obj.f.compute_prox(obj.A'*s/obj.mu, gam/obj.mu);
+        otherwise
+            error('not implemented');
+    end
+    p = obj.A*obj.x;
+end

+forbes/+functions/@HingeLoss/HingeLoss.m

@@ -0,0 +1,18 @@
+% HINGELOSS Hinge loss function
+
+classdef HingeLoss < forbes.functions.Proximable
+    properties
+        mu, b
+    end
+    methods
+        function obj = HingeLoss(mu, b)
+            if nargin < 1, mu = 1.0; end
+            if nargin < 2, b = 1.0; end
+            obj.mu = mu;
+            obj.b = b;
+        end
+        function p = is_convex(obj)
+            p = true;
+        end
+    end
+end

+forbes/+functions/@HingeLoss/compute_prox.m

@@ -0,0 +1,7 @@
+function [p, v] = compute_prox(obj, x, gam)
+    bx = obj.b .* x;
+    ind = bx < 1;
+    p(ind,1) = obj.b(ind) .* min(bx(ind)+gam*obj.mu,1);
+    p(~ind,1) = x(~ind);
+    v = obj.mu*sum(max(0,1-obj.b .* p));
+end

+forbes/+functions/@HuberLoss/HuberLoss.m

@@ -0,0 +1,18 @@
+% HUBERLOSS Huber loss function
+
+classdef HuberLoss < forbes.functions.Proximable
+    properties
+        del % coefficient
+    end
+    methods
+        function obj = HuberLoss(del)
+            obj.del = del;
+        end
+        function p = is_convex(obj)
+            p = true;
+        end
+        function p = is_smooth(obj)
+            p = true;
+        end
+    end
+end

+forbes/+functions/@HuberLoss/compute_gradient.m

@@ -0,0 +1,13 @@
+function [grad, val] = compute_gradient(obj, x)
+    absx = abs(x);
+    small = absx <= obj.del;
+    large = ~small;
+    sqx = (0.5/obj.del)*(x(small).^2);
+    linx = absx(large)-0.5*obj.del;
+    val = sum(sqx)+sum(linx);
+    if nargout >= 2
+        grad = zeros(length(x),1);
+        grad(small) = x(small)/obj.del;
+        grad(large) = sign(x(large));
+    end
+end

+forbes/+functions/@IndAffine/IndAffine.m

@@ -0,0 +1,15 @@
+% INDAFFINE Indicator function of an affine subspace
+
+classdef IndAffine < forbes.functions.Proximable
+    properties
+        A, b % define the affine subspace
+        L_prox % stores the Cholesky factor needed to compute prox
+    end
+    methods
+        function obj = IndAffine(A, b)
+            obj.A = A;
+            obj.b = b;
+            obj.L_prox = [];
+        end
+    end
+end

+forbes/+functions/@IndAffine/compute_prox.m

@@ -0,0 +1,8 @@
+function [p, v] = compute_prox(obj, x, gam)
+    if isempty(obj.L_prox)
+        obj.make_prox();
+    end
+    res = obj.A*x - obj.b;
+    p = x - obj.A'*(obj.L_prox'\(obj.L_prox\res));
+    v = 0.0;
+end

+forbes/+functions/@IndAffine/make_prox.m

@@ -0,0 +1,4 @@
+function make_prox(obj)
+    % can we make this more efficient, e.g., for sparse matrices?
+    obj.L_prox = chol(obj.A*obj.A','lower');
+end