SPAMS: a SPArse Modeling Software, v2.5

Contents

• 1  Introduction
• 2  Installation
• 3  Dictionary Learning and Matrix Factorization Toolbox
  • 3.1  Function mexTrainDL
  • 3.2  Function mexTrainDL_Memory
  • 3.3  Function mexStructTrainDL
  • 3.4  Function nmf
  • 3.5  Function nnsc
  • 3.6  Function mexArchetypalAnalysis
• 4  Sparse Decomposition Toolbox
  • 4.1  Function mexOMP
  • 4.2  Function mexOMPMask
  • 4.3  Function mexRidgeRegression
  • 4.4  Function mexLasso
  • 4.5  Function mexLassoWeighted
  • 4.6  Function mexLassoMask
  • 4.7  Function mexCD
  • 4.8  Function mexSOMP
  • 4.9  Function mexL1L2BCD
  • 4.10  Function mexSparseProject
  • 4.11  Function mexDecompSimplex
• 5  Proximal Toolbox
  • 5.1  Regularization Functions
  • 5.2  Function mexProximalFlat
  • 5.3  Function mexProximalTree
  • 5.4  Function mexProximalGraph
  • 5.5  Function mexProximalPathCoding
  • 5.6  Function mexEvalPathCoding
  • 5.7  Problems Addressed
    • 5.7.1  Regression Problems with the Square Loss
    • 5.7.2  Classification Problems with the Logistic Loss
    • 5.7.3  Multi-class Classification Problems with the Softmax Loss
    • 5.7.4  Multi-task Regression Problems with the Square Loss
    • 5.7.5  Multi-task Classification Problems with the Logistic Loss
    • 5.7.6  Multi-task and Multi-class Classification Problems with the Softmax Loss
  • 5.8  Function mexFistaFlat
  • 5.9  Function mexFistaTree
  • 5.10  Function mexFistaGraph
  • 5.11  Function mexFistaPathCoding
  • 5.12  Function solverPoisson
  • 5.13  Function mexIncrementalProx
  • 5.14  Function mexStochasticProx
• 6  A few Functions for manipulating images
  • 6.1  Function mexExtractPatches
  • 6.2  Function mexCombinePatches
  • 6.3  Function mexConvFistaFlat
• 7  Miscellaneous Functions
  • 7.1  Function mexConjGrad
  • 7.2  Function mexBayer
  • 7.3  Function mexCalcAAt
  • 7.4  Function mexCalcXAt
  • 7.5  Function mexCalcXY
  • 7.6  Function mexCalcXYt
  • 7.7  Function mexCalcXtY
  • 7.8  Function mexInvSym
  • 7.9  Function mexNormalize
  • 7.10  Function mexSort
  • 7.11  Function mexDisplayPatches
  • 7.12  Function mexCountPathsDAG
  • 7.13  Function mexRemoveCyclesGraph
  • 7.14  Function mexCountConnexComponents
  • 7.15  Function mexGraphOfGroupStruct
  • 7.16  Function mexGroupStructOfString
  • 7.17  Function mexReadGroupStruct
  • 7.18  Function mexTreeOfGroupStruct
  • 7.19  Function mexSimpleGroupTree
• A  Duality Gaps with Fenchel Duality
  • A.0.1  Duality Gaps without Intercepts
  • A.0.2  Duality Gaps with Intercepts

1  Introduction

SPAMS (SPArse Modeling Software) is an open-source optimization toolbox for sparse estimation with licence GPLv3. It implements algorithms for solving machine learning and signal processing problems involving sparse regularizations.

The library is coded in C++, is compatible with Linux, Mac, and Windows 32bits and 64bits Operating Systems. It is interfaced with Matlab, R and Python, but can be called from any C++ application (by hacking the code a bit).

It requires an implementation of BLAS and LAPACK for performing linear algebra operations. The ones shipped with matlab and R can be used, but also external libraries such as atlas, the netlib implementation, or the Intel Math Kernel Library can be used. It also exploits multi-core CPUs when this feature is supported by the compiler, through OpenMP.

The current licence is GPLv3, which is available at http://www.gnu.org/licenses/gpl.html. For other licensing possibilities allowing its use in proprietary softwares, please contact the author.

Version 2.5 of SPAMS is divided into several “toolboxes” and has a few additional miscellaneous functions:

• The Dictionary learning and matrix factorization toolbox contains the online learning technique of [20, 21] and its variants for solving various matrix factorization problems:
  • dictionary Learning for sparse coding;
  • sparse principal component analysis (seen as a sparse matrix factorization problem);
  • non-negative matrix factorization;
  • non-negative sparse coding;
  • dictionary learning with structured sparsity;
  • archetypal analysis [7, 37].
• The Sparse decomposition toolbox contains efficient implementations of
  • Orthogonal Matching Pursuit, (or Forward Selection) [35, 27];
  • the LARS/homotopy algorithm [30, 9] (variants for solving Lasso and Elastic-Net problems);
  • a weighted version of LARS;
  • OMP and LARS when data comes with a binary mask;
  • a coordinate-descent algorithm for ℓ₁-decomposition problems [12, 10, 36];
  • a greedy solver for simultaneous signal approximation as defined in [34, 33] (SOMP);
  • a solver for simulatneous signal approximation with ℓ₁/ℓ₂-regularization based on block-coordinate descent;
  • a homotopy method for the Fused-Lasso Signal Approximation as defined in [10] with the homotopy method presented in the appendix of [21];
  • a tool for projecting efficiently onto a few convex sets inducing sparsity such as the ℓ₁-ball using the method of [3, 18, 8], and Elastic-Net or Fused Lasso constraint sets as proposed in the appendix of [21].
  • an active-set algorithm for simplex decomposition problems [37].
• The Proximal toolbox: An implementation of proximal methods (ISTA and FISTA [1]) for solving a large class of sparse approximation problems with different combinations of loss and regularizations. One of the main features of this toolbox is to provide a robust stopping criterion based on duality gaps to control the quality of the optimization, whenever possible. It also handles sparse feature matrices for large-scale problems. The following regularizations are implemented:
  • Tikhonov regularization (squared ℓ₂-norm);
  • ℓ₁-norm, ℓ₂, ℓ_∞-norms;
  • Elastic-Net [39];
  • Fused Lasso [32];
  • tree-structured sum of ℓ₂-norms (see [15, 16]);
  • tree-structured sum of ℓ_∞-norms (see [15, 16]);
  • general sum of ℓ_∞-norms (see [22, 23]);
  • mixed ℓ₁/ℓ₂-norms on matrices [38, 29];
  • mixed ℓ₁/ℓ_∞-norms on matrices [38, 29];
  • mixed ℓ₁/ℓ₂-norms on matrices plus ℓ₁ [31, 11];
  • mixed ℓ₁/ℓ_∞-norms on matrices plus ℓ₁;
  • group-lasso with ℓ₂ or ℓ_∞-norms;
  • group-lasso+ℓ₁;
  • multi-task tree-structured sum of ℓ_∞-norms (see [22, 23]);
  • trace norm;
  • ℓ₀ pseudo-norm (only with ISTA);
  • tree-structured ℓ₀ (only with ISTA);
  • rank regularization for matrices (only with ISTA);
  • the path-coding penalties of [24].
  All of these regularization functions can be used with the following losses
  • square loss;
  • square loss with missing observations;
  • logistic loss, weighted logistic loss;
  • multi-class logistic.
  This toolbox can also enforce non-negativity constraints, handle intercepts and sparse matrices. There are also a few additional undocumented functionalities, which are available in the source code. For some combinations of loss and regularizers, stochastic and incremental proximal gradient solvers are also implemented [26, 25].
• A few tools for performing linear algebra operations such as a conjugate gradient algorithm, manipulating sparse matrices and graphs.

The toolbox was written by Julien Mairal at INRIA, with the collaboration of Francis Bach (INRIA), Jean Ponce (Ecole Normale Supérieure), Guillermo Sapiro (University of Minnesota), Guillaume Obozinski (INRIA) and Rodolphe Jenatton (INRIA).

R and Python interfaces have been written by Jean-Paul Chieze (INRIA). The archetypal analysis implementation was written by Yuansi Chen, during an internship at INRIA, with the collaboration of Zaid Harchaoui.

2  Installation

The toolbox used to come with pre-compiled binaries for various plateforms. Now the sources are available and you will have to compile it yourself.

The user has the choice of the BLAS library, but the Intel MKL is recommended for the best performance. Note that the builtin blas library of Matlab is a version of the Intel MKL (not the most recent one though), and offers a good performance.

The folder doc contains the documentation in pdf and html. build/ contains the binary files, including the help for each command. test_release contains various matlab scripts which call the different functions of this toolbox.

The software package comes also with a script bash that has to be used to launch matlab for Linux and/or Mac OS versions. The install procedure is described in a file called HOW_TO_INSTALL and in the file compile.m for the Matlab version.

3  Dictionary Learning and Matrix Factorization Toolbox

This is the section for dictionary learning and matrix factorization, corresponding to [20, 21].

3.1  Function mexTrainDL

This is the main function of the toolbox, implementing the learning algorithms of [21]. Given a training set x¹,…, . It aims at solving

ψ is a sparsity-inducing regularizer and C is a constraint set for the dictionary. As shown in [21] and in the help file below, various combinations can be used for ψ and C for solving different matrix factorization problems. What is more, positivity constraints can be added to α as well. The function admits several modes for choosing the optimization parameters, using the parameter-free strategy proposed in [20], or using the parameters t₀ and ρ presented in [21]. Note that for problems of a reasonable size, and when ψ is the ℓ₁-norm, the function mexTrainDL_Memory can be faster but uses more memory.

The following piece of code illustrates how to use this function.

3.2  Function mexTrainDL_Memory

Memory-consuming version of mexTrainDL. This function is well adapted to small/medium-size problems: It requires storing all the coefficients α and is therefore impractical for very large datasets. However, in many situations, one can afford this memory cost and it is better to use this method, which is faster than mexTrainDL. Note that unlike mexTrainDL this function does not allow warm-restart.

The following piece of code illustrates how to use this function.

3.3  Function mexStructTrainDL

This function allows to use mexTrainDL with structured regularization functions for the coefficients α. It internally uses the FISTA algorithm.

The following piece of code illustrates how to use this function.

3.4  Function nmf

This function is an example on how to use the function mexTrainDL for the problem of non-negative matrix factorization formulated in [17]. Note that mexTrainDL can be replaced by mexTrainDL_Memory in this function for small or medium datasets.

The following piece of code illustrates how to use this function.

3.5  Function nnsc

This function is an example on how to use the function mexTrainDL for the problem of non-negative sparse coding as defined in [14]. Note that mexTrainDL can be replaced by mexTrainDL_Memory in this function for small or medium datasets.

3.6  Function mexArchetypalAnalysis

The function optimizes the archetypal analysis formulation of [7]. It follows the methodology presented in the following paper [37].

The following piece of code illustrates how to use this function.

4  Sparse Decomposition Toolbox

This toolbox implements several algorithms for solving signal reconstruction problems. It is mostly adapted for solving a large number of small/medium scale problems, but can be also efficient sometimes with large scale ones.

4.1  Function mexOMP

This is a fast implementation of the Orthogonal Matching Pursuit algorithm (or forward selection) [27, 35]. Given a matrix of signals X=[x¹,…,xⁿ] in ℝ^{m × n} and a dictionary D=[d¹,…,d^p] in ℝ^{m × p}, the algorithm computes a matrix A=[α¹,…,αⁿ] in ℝ^{p × n}, where for each column x of X, it returns a coefficient vector α which is an approximate solution of the following NP-hard problem

or

or

For efficienty reasons, the method first computes the covariance matrix D^TD, then for each signal, it computes D^Tx and performs the decomposition with a Cholesky-based algorithm (see [6] for instance).

Note that mexOMP can return the “greedy” regularization path if needed (see below):

The following piece of code illustrates how to use this function.

4.2  Function mexOMPMask

This is a variant of mexOMP with the possibility of handling a binary mask. Given a binary mask B=[β¹,…,βⁿ] in {0,1}^{m × n}, it returns a matrix A=[α¹,…,αⁿ] such that for every column x of X, β of B, it computes a column α of A by addressing

or

or

where diag(β) is a diagonal matrix with the entries of β on the diagonal.

The following piece of code illustrates how to use this function.

4.3  Function mexRidgeRegression

This is a ridge regression solver using a conjugate gradient solver.

The following piece of code illustrates how to use this function.

4.4  Function mexLasso

This is a fast implementation of the LARS algorithm [9] (variant for solving the Lasso) for solving the Lasso or Elastic-Net. Given a matrix of signals X=[x¹,…,xⁿ] in ℝ^{m × n} and a dictionary D in ℝ^{m × p}, depending on the input parameters, the algorithm returns a matrix of coefficients A=[α¹,…,αⁿ] in ℝ^{p × n} such that for every column x of X, the corresponding column α of A is the solution of

or

or

For efficiency reasons, the method first compute the covariance matrix D^TD, then for each signal, it computes D^Tx and performs the decomposition with a Cholesky-based algorithm (see [9] for instance). The implementation has also an option to add positivity constraints on the solutions α. When the solution is very sparse and the problem size is reasonable, this approach can be very efficient. Moreover, it gives the solution with an exact precision, and its performance does not depend on the correlation of the dictionary elements, except when the solution is not unique (the algorithm breaks in this case).

Note that mexLasso can return the whole regularization path of the first signal x₁ and can handle implicitely the matrix D if the quantities D^TD and D^Tx are passed as an argument, see below:

The following piece of code illustrates how to use this function.

4.5  Function mexLassoWeighted

This is a fast implementation of a weighted version of LARS [9]. Given a matrix of signals X=[x¹,…,xⁿ] in ℝ^{m × n}, a matrix of weights W=[w¹,…,wⁿ] ∈ ℝ^{p × n}, and a dictionary D in ℝ^{m × p}, depending on the input parameters, the algorithm returns a matrix of coefficients A=[α¹,…,αⁿ] in ℝ^{p × n}, such that for every column x of X, w of W, it computes a column α of A, which is the solution of

or

or

The implementation has also an option to add positivity constraints on the solutions α. This function is potentially useful for implementing efficiently the randomized Lasso of [28], or reweighted-ℓ₁ schemes [4].

The following piece of code illustrates how to use this function.

4.6  Function mexLassoMask

This is a variant of mexLasso with the possibility of adding a mask B=[β¹,…,βⁿ], as in mexOMPMask. For every column x of X, β of B, it computes a column α of A, which is the solution of

or

or

The following piece of code illustrates how to use this function.

4.7  Function mexCD

Coordinate-descent approach for solving Eq. (10) and Eq. (9). Note that unlike mexLasso, it is not implemented to solve the Elastic-Net formulation. To solve Eq. (9), the algorithm solves a sequence of problems of the form (10) using simple heuristics. Coordinate descent is very simple and in practice very powerful. It performs better when the correlation between the dictionary elements is small.

The following piece of code illustrates how to use this function.

4.8  Function mexSOMP

This is a fast implementation of the Simultaneous Orthogonal Matching Pursuit algorithm. Given a set of matrices X=[X¹,…,Xⁿ] in ℝ^{m × N}, where the Xⁱ’s are in ℝ^{m × n_i}, and a dictionary D in ℝ^{m × p}, the algorithm returns a matrix of coefficients A=[A¹,…,Aⁿ] in ℝ^{p × N} which is an approximate solution of the following NP-hard problem

or

To be efficient, the method first compute the covariance matrix D^TD, then for each signal, it computes D^TXⁱ and performs the decomposition with a Cholesky-based algorithm.

The following piece of code illustrates how to use this function.

4.9  Function mexL1L2BCD

This is a fast implementation of a simultaneous signal decomposition formulation. Given a set of matrices X=[X¹,…,Xⁿ] in ℝ^{m × N}, where the Xⁱ’s are in ℝ^{m × n_i}, and a dictionary D in ℝ^{m × p}, the algorithm returns a matrix of coefficients A=[A¹,…,Aⁿ] in ℝ^{p × N} which is an approximate solution of the following NP-hard problem

or

To be efficient, the method first compute the covariance matrix D^TD, then for each signal, it computes D^TXⁱ and performs the decomposition with a Cholesky-based algorithm.

The following piece of code illustrates how to use this function.

4.10  Function mexSparseProject

This is a multi-purpose function, implementing fast algorithms for projecting on convex sets, but it also solves the fused lasso signal approximation problem. The proposed method is detailed in [21]. The main problems addressed by this function are the following: Given a matrix U=[u₁,…,u_n] in ℝ^{m × n}, it finds a matrix V=[v₁,…,v_n] in ℝ^{m × n} so that for all column u of U, it computes a column v of V solving

or

or

or

Note that for the two last cases, the method performs a small approximation. The method follows the regularization path, goes from one kink to another, and stop whenever the constraint is not satisfied anymore. The solution returned by the algorithm is the one obtained at the last kink of the regularization path, which is in practice close, but not exactly the same as the solution. This will be corrected in a future release of the toolbox.

The following piece of code illustrates how to use this function.

4.11  Function mexDecompSimplex

The function implements an active-set algorithm [37] for solving

5  Proximal Toolbox

The previous toolbox we have presented is well adapted for solving a large number of small and medium-scale sparse decomposition problems with the square loss, which is typical from the classical dictionary learning framework. We now present a new software package that is adapted for solving a wide range of possibly large-scale learning problems, with several combinations of losses and regularization terms. The method implements the proximal methods of [1], and includes the proximal solvers for the tree-structured regularization of [15], and the solver of [22] for general structured sparse regularization. The solver for structured sparse regularization norms includes a C++ max-flow implementation of the push-relabel algorithm of [13], with heuristics proposed by [5].

This implementation also provides robust stopping criteria based on duality gaps, which are presented in Appendix A. It can handle intercepts (unregularized variables). The general formulation that our software can solve take the form

where f is a smooth loss function and ψ is a regularization function. When one optimizes a matrix W in ℝ^{p × r} instead of a vector w in ℝ^p, we will write

Note that the software can possibly handle nonnegativity constraints.

We start by presenting the type of regularization implemented in the software

5.1  Regularization Functions

Our software can handle the following regularization functions ψ for vectors w in ℝ^p:

• The Tikhonov regularization: ψ(w) =^▵ 1/2||w||₂².
• The ℓ₁-norm: ψ(w) =^▵ ||w||₁.
• The Elastic-Net: ψ(w) =^▵ ||w||₁+γ||w||₂².
• The Fused-Lasso: ψ(w) =^▵ ||w||₁+γ||w||₂²+γ₂∑_i=1^p−1|w_i+1−w_i|.
• The group Lasso: ψ(w) =^▵ ∑_{g ∈ G} η_g ||w_g||₂, where G are groups of variables.
• The group Lasso with ℓ_∞-norm: ψ(w) =^▵ ∑_{g ∈ G} η_g ||w_g||_∞, where G are groups of variables.
• The sparse group Lasso: same as above but with an additional ℓ₁ term.
• The tree-structured sum of ℓ₂-norms: ψ(w) =^▵ ∑_{g ∈ G} η_g ||w_g||₂, where G is a tree-structured set of groups [15], and the η_g are positive weights.
• The tree-structured sum of ℓ_∞-norms: ψ(w) =^▵ ∑_{g ∈ G} η_g ||w_g||_∞. See [15]
• General sum of ℓ_∞-norms: ψ(w) =^▵ ∑_{g ∈ G} η_g ||w_g||_∞, where no assumption are made on the groups G.
• The path-coding penalties of [24].
• the ℓ₁-constraint.
• a few undocumented ones.

Our software also handles regularization functions ψ on matrices W in ℝ^{p × r} (note that W can be transposed in these formulations). In particular,

• The ℓ₁/ℓ₂-norm: ψ(W) =^▵ ∑_i=1^p ||W_i||₂, where W_i denotes the i-th row of W.
• The ℓ₁/ℓ_∞-norm: ψ(W) =^▵ ∑_i=1^p ||W_i||_∞,
• The ℓ₁/ℓ₂+ℓ₁-norm: ψ(W) =^▵ ∑_i=1^p ||W_i||₂ + λ₂ ∑_i,j|W_ij|.
• The ℓ₁/ℓ_∞+ℓ₁-norm: ψ(W) =^▵ ∑_i=1^p ||W_i||_∞+λ₂ ∑_i,j|W_ij|,
• The ℓ₁/ℓ_∞-norm on rows and columns: ψ(W) =^▵ ∑_i=1^p ||W_i||_∞+λ₂ ∑_j=1^r||W^j||_∞, where W^j denotes the j-th column of W.
• The multi-task tree-structured sum of ℓ_∞-norms:

  ψ(W)

  ▵ =

  r ∑ i=1

  ∑ g ∈ G

  ηg||wgi||∞+ γ

  ∑ g ∈ G

  ηg

  max j ∈ g

  ||Wj||∞,     (25)

  where the first double sums is in fact a sum of independent structured norms on the columns wⁱ of W, and the right term is a tree-structured regularization norm applied to the ℓ_∞-norm of the rows of W, thereby inducing the tree-structured regularization at the row level. G is here a tree-structured set of groups.
• The multi-task general sum of ℓ_∞-norms is the same as Eq. (25) except that the groups G are general overlapping groups.
• The trace norm: ψ(W) =^▵ ||W||_*.

Non-convex regularizations are also implemented with the ISTA algorithm (no duality gaps are of course provided in these cases):

• The ℓ₀-pseudo-norm: ψ(w) =^▵ ||w||₀.
• The rank: ψ(W) =^▵ randk(W).
• The tree-structured ℓ₀-pseudo-norm: ψ(w) =^▵ ∑_{g ∈ G} δ_{wg ≠ 0}.

All of these regularization terms for vectors or matrices can be coupled with nonnegativity constraints. It is also possible to add an intercept, which one wishes not to regularize, and we will include this possibility in the next sections. There are also a few hidden undocumented options which are available in the source code.

We now present 3 functions for computing proximal operators associated to the previous regularization functions.

5.2  Function mexProximalFlat

This function computes the proximal operators associated to many regularization functions, for input signals U=[u¹,…,uⁿ] in ℝ^{p × n}, it finds a matrix V=[v¹,…,vⁿ] in ℝ^{p × n} such that:

•  If one chooses a regularization function on vectors, for every column u of U, it computes one column v of V solving

or

or

or

or

or

where T is a tree-structured set of groups (see [16]), and δ^g(v) = 0 if v_g=0 and 1 otherwise. It can also solve

or

or

where G is any kind of set of groups.

This function can also solve the following proximal operators on matrices

where V_i is the i-th row of V, or

or

or

or

where V^j is the j-th column of V.

See details below:

The following piece of code illustrates how to use this function.

5.3  Function mexProximalTree

This function computes the proximal operators associated to tree-structured regularization functions, for input signals U=[u¹,…,uⁿ] in ℝ^{p × n}, and a tree-structured set of groups [15], it computes a matrix V=[v¹,…,vⁿ] in ℝ^{p × n}. When one uses a regularization function on vectors, it computes a column v of V for every column u of U:

or

or

where δ^g(v)=0 if v_g=0 and 1 otherwise (see appendix of [16]).

When the multi-task tree-structured regularization function is used, it solves

which is a formulation presented in [22].

This function can also be used for computing the proximal operators addressed by mexProximalFlat (it will just not take into account the tree structure). The way the tree is incoded is presented below, (and examples are given in the file test_ProximalTree.m, with more usage details:

The following piece of code illustrates how to use this function.

5.4  Function mexProximalGraph

This function computes the proximal operators associated to structured sparse regularization, for input signals U=[u¹,…,uⁿ] in ℝ^{p × n}, and a set of groups [22], it returns a matrix V=[v¹,…,vⁿ] in ℝ^{p × n}. When one uses a regularization function on vectors, it computes a column v of V for every column u of U:

or with a regularization function on matrices, it computes V solving

This function can also be used for computing the proximal operators addressed by mexProximalFlat. The way the graph is incoded is presented below (and also in the example file test_ProximalGraph.m, with more usage details:

The following piece of code illustrates how to use this function.

5.5  Function mexProximalPathCoding

This function computes the proximal operators associated to the path coding penalties of [24].

The following piece of code illustrates how to use this function.

This function is associated to a function to evaluate the penalties:

5.6  Function mexEvalPathCoding

The following piece of code illustrates how to use this function.

After having presented the regularization terms which our software can handle, we present the various formulations that we address

5.7  Problems Addressed

We present here regression or classification formulations and their multi-task variants.

5.7.1  Regression Problems with the Square Loss

Given a training set {xⁱ,y_i}_i=1ⁿ, with xⁱ ∈ ℝ^p and y_i ∈ ℝ for all i in [ 1;n ], we address

where b is an optional variable acting as an “intercept”, which is not regularized, and ψ can be any of the regularization functions presented above. Let us consider the vector y in ℝⁿ that carries the entries y_i. The problem without the intercept takes the following form, which we have already encountered in the previous toolbox, but with different notations:

where the X=[xⁱ,…,xⁿ]^T (the xⁱ’s are here the rows of X).

5.7.2  Classification Problems with the Logistic Loss

The next formulation that our software can solve is the regularized logistic regression formulation. We are again given a training set {xⁱ,y_i}_i=1ⁿ, with xⁱ ∈ ℝ^p, but the variables y_i are now in {−1,+1} for all i in [ 1;n ]. The optimization problem we address is

with again ψ taken to be one of the regularization function presented above, and b is an optional intercept.

5.7.3  Multi-class Classification Problems with the Softmax Loss

We have also implemented a multi-class logistic classifier (or softmax). For a classification problem with r classes, we are given a training set {xⁱ,y_i}_i=1ⁿ, where the variables xⁱ are still vectors in ℝ^p, but the y_i’s have integer values in {1,2,…,r}. The formulation we address is the following multi-class learning problem

where W = [w¹,…,w^r] and the optional vector b in ℝ^r carries intercepts for each class.

5.7.4  Multi-task Regression Problems with the Square Loss

We are now considering a problem with r tasks, and a training set {xⁱ,yⁱ}_i=1ⁿ, where the variables xⁱ are still vectors in ℝ^p, and yⁱ is a vector in ℝ^r. We are looking for r regression vectors w^j, for j∈ [ 1;r ], or equivalently for a matrix W=[w¹,…,w^r] in ℝ^{p × r}. The formulation we address is the following multi-task regression problem

where ψ is any of the regularization function on matrices we have presented in the previous section. Note that by introducing the appropriate variables Y, the problem without intercept could be equivalently rewritten

5.7.5  Multi-task Classification Problems with the Logistic Loss

The multi-task version of the logistic regression follows the same principle. We consider r tasks, and a training set {xⁱ,yⁱ}_i=1ⁿ, with the xⁱ’s in ℝ^p, and the yⁱ’s are vectors in {−1,+1}^r. We look for a matrix W=[w¹,…,w^r] in ℝ^{p × r}. The formulation is the following multi-task regression problem

5.7.6  Multi-task and Multi-class Classification Problems with the Softmax Loss

The multi-task/multi-class version directly follows from the formulation of Eq. (46), but associates with each class a task, and as a consequence, regularizes the matrix W in a particular way:

How duality gaps are computed for any of these formulations is presented in Appendix A. We now present the main functions for solving these problems

5.8  Function mexFistaFlat

Given a matrix X=[x¹,…,x^p]^T in ℝ^{m × p}, and a matrix Y=[y¹,…,yⁿ], it solves the optimization problems presented in the previous section, with the same regularization functions as mexProximalFlat. see usage details below:

The following piece of code illustrates how to use this function.

5.9  Function mexFistaTree

Given a matrix X=[x¹,…,x^p]^T in ℝ^{m × p}, and a matrix Y=[y¹,…,yⁿ], it solves the optimization problems presented in the previous section, with the same regularization functions as mexProximalTree. see usage details below: