Previous Up Next

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

 
min
w ∈ ℝp
 [g(w
=
 
 f(w) + λψ(w)],

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

 
min
W ∈ ℝp × r
 [g(W
=
 
 f(W) + λψ(W)].

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:

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

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

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 spams.proximalFlat

This function computes the proximal operators associated to many regularization functions, for input signals U=[u1,…,un] in ℝp × n, it finds a matrix V=[v1,…,vn] 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

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ ||v||0,     (26)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ ||v||1,     (27)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + 
λ
2
 ||v||22,     (28)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ ||v||1 + λ2||v||22,     (29)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ
p−1
j=1
|vj+1ivji|+λ2 ||v||1 + λ3||v||22,     (30)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ T
 δg(v),     (31)

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

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ T
 ηg ||vg||2,     (32)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ T
 ηg ||vg||,     (33)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ G
 ηg ||vg||,     (34)

where G is any kind of set of groups.

This function can also solve the following proximal operators on matrices

 
min
V ∈ ℝp × n
 
1
2
||UV||F2 + λ 
p
i=1
 ||Vi||2,     (35)

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

 
min
V ∈ ℝp × n
 
1
2
||UV||F2 + λ 
p
i=1
 ||Vi||,     (36)

or

 
min
V ∈ ℝp × n
 
1
2
||UV||F2 + λ 
p
i=1
 ||Vi||2 +λ2 
p
i=1
n
j=1
 |Vij|,     (37)

or

 
min
V ∈ ℝp × n
 
1
2
||UV||F2 + λ 
p
i=1
 ||Vi||2 
p
i=1
n
j=1
 |Vij|,     (38)

or

 
min
V ∈ ℝp × n
 
1
2
||UV||F2 + λ 
p
i=1
 ||Vi||2 
n
j=1
 ||Vj||.     (39)

where Vj is the j-th column of V.

See details below:

#
# Name: spams.proximalFlat
#
# Usage: spams.proximalFlat(U,return_val_loss = FALSE,numThreads =-1,lambda1=1.0,lambda2=0.,
#                    lambda3=0.,intercept=FALSE,resetflow=FALSE,regul="",verbose=FALSE,
#                    pos=FALSE,clever=TRUE,size_group=1,groups = NULL,transpose=FALSE)
#
# Description:
#     spams.proximalFlat computes proximal operators. Depending
#         on the value of regul, it computes 
#         
#         Given an input matrix U=[u^1,\ldots,u^n], it computes a matrix 
#         V=[v^1,\ldots,v^n] such that
#         if one chooses a regularization functions on vectors, it computes
#         for each column u of U, a column v of V solving
#         if regul='l0'
#             argmin 0.5||u-v||_2^2 + lambda1||v||_0
#         if regul='l1'
#             argmin 0.5||u-v||_2^2 + lambda1||v||_1
#         if regul='l2'
#             argmin 0.5||u-v||_2^2 + 0.5lambda1||v||_2^2
#         if regul='elastic-net'
#             argmin 0.5||u-v||_2^2 + lambda1||v||_1 + lambda1_2||v||_2^2
#         if regul='fused-lasso'
#             argmin 0.5||u-v||_2^2 + lambda1 FL(v) + ...
#                               ...  lambda1_2||v||_1 + lambda1_3||v||_2^2
#         if regul='linf'
#             argmin 0.5||u-v||_2^2 + lambda1||v||_inf
#         if regul='l1-constraint'
#             argmin 0.5||u-v||_2^2 s.t. ||v||_1 <= lambda1
#         if regul='l2-not-squared'
#             argmin 0.5||u-v||_2^2 + lambda1||v||_2
#         if regul='group-lasso-l2'  
#             argmin 0.5||u-v||_2^2 + lambda1 sum_g ||v_g||_2 
#             where the groups are either defined by groups or by size_group,
#         if regul='group-lasso-linf'
#             argmin 0.5||u-v||_2^2 + lambda1 sum_g ||v_g||_inf
#         if regul='sparse-group-lasso-l2'  
#             argmin 0.5||u-v||_2^2 + lambda1 sum_g ||v_g||_2 + lambda1_2 ||v||_1
#             where the groups are either defined by groups or by size_group,
#         if regul='sparse-group-lasso-linf'
#             argmin 0.5||u-v||_2^2 + lambda1 sum_g ||v_g||_inf + lambda1_2 ||v||_1
#         if regul='trace-norm-vec' 
#             argmin 0.5||u-v||_2^2 + lambda1 ||mat(v)||_* 
#            where mat(v) has size_group rows
#            
#         if one chooses a regularization function on matrices
#         if regul='l1l2',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1||V||_{1/2}
#         if regul='l1linf',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1||V||_{1/inf}
#         if regul='l1l2+l1',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1||V||_{1/2} + lambda1_2||V||_{1/1}
#         if regul='l1linf+l1',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1||V||_{1/inf} + lambda1_2||V||_{1/1}
#         if regul='l1linf+row-column',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1||V||_{1/inf} + lambda1_2||V'||_{1/inf}
#         if regul='trace-norm',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1||V||_*
#         if regul='rank',  V= 
#             argmin 0.5||U-V||_F^2 + lambda1 rank(V)
#         if regul='none',  V= 
#             argmin 0.5||U-V||_F^2 
#             
#         for all these regularizations, it is possible to enforce non-negativity constraints
#         with the option pos, and to prevent the last row of U to be regularized, with
#         the option intercept
#
# Inputs:
#       U:  double m x n matrix   (input signals)
#               m is the signal size
#       return_val_loss:     
#               if true the function will return 2 matrices in a list.
#       lambda1:  (regularization parameter)
#       regul: (choice of regularization, see above)
#       lambda2:  (optional, regularization parameter)
#       lambda3:  (optional, regularization parameter)
#       verbose: (optional, verbosity level, false by default)
#       intercept: (optional, last row of U is not regularized,
#         false by default)
#       transpose: (optional, transpose the matrix in the regularization function)
#       size_group: (optional, for regularization functions assuming a group
#         structure). It is a scalar. When groups is not specified, it assumes
#         that the groups are the sets of consecutive elements of size size_group
#       groups: (int32, optional, for regularization functions assuming a group
#         structure. It is an int32 vector of size m containing the group indices of the
#         variables (first group is 1).
#       pos: (optional, adds positivity constraints on the
#         coefficients, false by default)
#       numThreads: (optional, number of threads for exploiting
#         multi-core / multi-cpus. By default, it takes the value -1,
#         which automatically selects all the available CPUs/cores).
#       resetflow:    undocumented; modify at your own risks!
#       clever:    undocumented; modify at your own risks!
#
# Output:
#       V: double m x n matrix (output coefficients)
#       val_regularizer: double 1 x n vector (value of the regularization
#       term at the optimum).
#       val_loss:        vector of size ncol(U)
#             alpha <- spams.proximalFlat(U,return_val_loss = FALSE, ...)
#             v <- spams.proximalFlat(U,return_val_loss = TRUE, ...)
#             alpha <- v[[1]]
#             val_loss <- v[[2]]
#
# Authors:
# Julien MAIRAL, 2010 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (R interface)
#
# Note:
#     Valid values for the regularization parameter (regul) are:
#       "l0", "l1", "l2", "linf", "l2-not-squared", "elastic-net", "fused-lasso",
#       "group-lasso-l2", "group-lasso-linf", "sparse-group-lasso-l2",
#       "sparse-group-lasso-linf", "l1l2", "l1linf", "l1l2+l1", "l1linf+l1",
#       "tree-l0", "tree-l2", "tree-linf", "graph", "graph-ridge", "graph-l2",
#       "multi-task-tree", "multi-task-graph", "l1linf-row-column", "trace-norm",
#       "trace-norm-vec", "rank", "rank-vec", "none"
#

The following piece of code illustrates how to use this function. The following piece of code contains usage examples:

library(spams)
m = 100;n = 1000
U = matrix(rnorm(m * n),nrow = m,ncol = n,byrow = FALSE)

# test l0
.printf("\nprox l0\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l0', pos = FALSE,intercept = FALSE)),n = 1)

# test l1
.printf("\nprox l1, intercept, positivity constraint\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1', pos = TRUE,intercept = TRUE)),n = 1)

# test l2
.printf("\nprox squared-l2\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l2', pos = FALSE,intercept = FALSE)),n = 1)

# test elastic-net
.printf("\nprox elastic-net\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'elastic-net', pos = FALSE,intercept = FALSE,lambda2 = 0.1)),n = 1)

# test fused-lasso
.printf("\nprox fused-lasso\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'fused-lasso', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

# test l1l2
.printf("\nprox mixed norm l1/l2\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1l2', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

# test l1linf
.printf("\nprox mixed norm l1/linf\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1linf', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

# test l1l2+l1
.printf("\nprox mixed norm l1/l2 + l1\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1l2+l1', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

# test l1linf+l1
.printf("\nprox mixed norm l1/linf + l1\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1linf+l1', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

# test l1linf-row-column
.printf("\nprox mixed norm l1/linf on rows and columns\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1linf-row-column', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

# test none
.printf("\nprox no regularization\n")
alpha = Xtest1('spams',quote(spams.proximalFlat(U,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'none', pos = FALSE,intercept = FALSE,lambda2 = 0.1,lambda3 = 0.1)),n = 1)

5.3  Function spams.proximalTree

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

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ T
 ηg ||vg||2,     (40)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ T
 ηg ||vg||,     (41)

or

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ T
 δg(v),     (42)

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

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

 
min
V ∈ ℝp× n
 
1
2
||UV||F2 + λ 
n
i=1
 
g ∈ T
 ηg ||vgi||+ λ2
 
g ∈ T
 
 
max
j ∈ g
 ||vgj||,     (43)

which is a formulation presented in [22].

This function can also be used for computing the proximal operators addressed by spams.proximalFlat (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:

#
# Name: spams.proximalTree
#
# Usage: spams.proximalTree(U,tree,return_val_loss = FALSE,numThreads =-1,lambda1=1.0,lambda2=0.,
#                    lambda3=0.,intercept=FALSE,resetflow=FALSE,regul="",verbose=FALSE,
#                    pos=FALSE,clever=TRUE,size_group=1,transpose=FALSE)
#
# Description:
#     spams.proximalTree computes a proximal operator. Depending
#         on the value of regul, it computes 
#         
#         Given an input matrix U=[u^1,\ldots,u^n], and a tree-structured set of groups T,
#         it returns a matrix V=[v^1,\ldots,v^n]:
#         
#         when the regularization function is for vectors,
#         for every column u of U, it compute a column v of V solving
#         if regul='tree-l0'
#             argmin 0.5||u-v||_2^2 + lambda1 \sum_{g \in T} \delta^g(v)
#         if regul='tree-l2'
#           for all i, v^i = 
#             argmin 0.5||u-v||_2^2 + lambda1\sum_{g \in T} \eta_g||v_g||_2
#         if regul='tree-linf'
#           for all i, v^i = 
#             argmin 0.5||u-v||_2^2 + lambda1\sum_{g \in T} \eta_g||v_g||_inf
#             
#         when the regularization function is for matrices:
#         if regul='multi-task-tree'
#            V=argmin 0.5||U-V||_F^2 + lambda1 \sum_{i=1}^n\sum_{g \in T} \eta_g||v^i_g||_inf + ...
#                                                lambda1_2 \sum_{g \in T} \eta_g max_{j in g}||V_j||_{inf}
#                                                
#         it can also be used with any non-tree-structured regularization addressed by spams.proximalFlat
#         
#         for all these regularizations, it is possible to enforce non-negativity constraints
#         with the option pos, and to prevent the last row of U to be regularized, with
#         the option intercept
#
# Inputs:
#       U:  double m x n matrix   (input signals)
#               m is the signal size
#       tree: named list 
#             with four fields, eta_g, groups, own_variables and N_own_variables.
#             
#             The tree structure requires a particular organization of groups and variables
#                * Let us denote by N = |T|, the number of groups.
#                  the groups should be ordered T={g1,g2,\ldots,gN} such that if gi is included
#                  in gj, then j <= i. g1 should be the group at the root of the tree 
#                  and contains every variable.
#                * Every group is a set of  contiguous indices for instance 
#                  gi={3,4,5} or gi={4,5,6,7} or gi={4}, but not {3,5};
#                * We define root(gi) as the indices of the variables that are in gi,
#                  but not in its descendants. For instance for
#                  T={ g1={1,2,3,4},g2={2,3},g3={4} }, then, root(g1)={1}, 
#                  root(g2)={2,3}, root(g3)={4},
#                  We assume that for all i, root(gi) is a set of contigous variables
#                * We assume that the smallest of root(gi) is also the smallest index of gi.
#                
#                For instance, 
#                  T={ g1={1,2,3,4},g2={2,3},g3={4} }, is a valid set of groups.
#                  but we can not have
#                  T={ g1={1,2,3,4},g2={1,2},g3={3} }, since root(g1)={4} and 4 is not the
#                  smallest element in g1.
#                  
#             We do not lose generality with these assumptions since they can be fullfilled for any
#             tree-structured set of groups after a permutation of variables and a correct ordering of the
#             groups.
#             see more examples in test_ProximalTree.m of valid tree-structured sets of groups.
#             
#             The first fields sets the weights for every group
#                tree[['eta_g']]            double N vector 
#                
#             The next field sets inclusion relations between groups 
#             (but not between groups and variables):
#                tree[['groups']]           sparse (double or boolean) N x N matrix  
#                the (i,j) entry is non-zero if and only if i is different than j and 
#                gi is included in gj.
#                the first column corresponds to the group at the root of the tree.
#                
#             The next field define the smallest index of each group gi, 
#             which is also the smallest index of root(gi)
#             tree[['own_variables']]    int32 N vector
#             
#             The next field define for each group gi, the size of root(gi)
#             tree[['N_own_variables']]  int32 N vector 
#             
#             examples are given in test_ProximalTree.m
#             
#       return_val_loss:     
#               if true the function will return 2 matrices in a list.
#       lambda1:  (regularization parameter)
#       regul: (choice of regularization, see above)
#       lambda2:  (optional, regularization parameter)
#       lambda3:  (optional, regularization parameter)
#       verbose: (optional, verbosity level, false by default)
#       intercept: (optional, last row of U is not regularized,
#         false by default)
#       pos: (optional, adds positivity constraints on the
#         coefficients, false by default)
#       transpose: (optional, transpose the matrix in the regularization function)
#       size_group: (optional, for regularization functions assuming a group
#         structure). It is a scalar. When groups is not specified, it assumes
#         that the groups are the sets of consecutive elements of size size_group
#       numThreads: (optional, number of threads for exploiting
#         multi-core / multi-cpus. By default, it takes the value -1,
#         which automatically selects all the available CPUs/cores).
#       resetflow:    undocumented; modify at your own risks!
#       clever:    undocumented; modify at your own risks!
#
# Output:
#       V: double m x n matrix (output coefficients)
#       val_regularizer: double 1 x n vector (value of the regularization
#       term at the optimum).
#       val_loss:        vector of size ncol(U)
#             alpha <- spams.proximalTree(U,tree,return_val_loss = FALSE, ...)
#             v <- spams.proximalTree(U,tree,return_val_loss = TRUE, ...)
#             alpha <- v[[1]]
#             val_loss <- v[[2]]
#
# Authors:
# Julien MAIRAL, 2010 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (R interface)
#
# Note:
#     Valid values for the regularization parameter (regul) are:
#       "l0", "l1", "l2", "linf", "l2-not-squared", "elastic-net", "fused-lasso",
#       "group-lasso-l2", "group-lasso-linf", "sparse-group-lasso-l2",
#       "sparse-group-lasso-linf", "l1l2", "l1linf", "l1l2+l1", "l1linf+l1",
#       "tree-l0", "tree-l2", "tree-linf", "graph", "graph-ridge", "graph-l2",
#       "multi-task-tree", "multi-task-graph", "l1linf-row-column", "trace-norm",
#       "trace-norm-vec", "rank", "rank-vec", "none"
#

The following piece of code illustrates how to use this function. The following piece of code contains usage examples:

library(spams)
m = 10;n = 1000
U = matrix(rnorm(m * n),nrow = m,ncol = n,byrow = FALSE)

.printf("First tree example\n")
# Example 1 of tree structure
# tree structured groups:
# g1= {0 1 2 3 4 5 6 7 8 9}
# g2= {2 3 4}
# g3= {5 6 7 8 9}

own_variables = as.vector(c(0,2,5),mode= 'integer')
N_own_variables = as.vector(c(2,3,5),mode= 'integer')
eta_g = as.vector(c(1,1,1),mode = 'double')
groups = matrix(as.vector(c(0,0,0,
  1,0,0,
  1,0,0),mode='logical'),ncol = 3,byrow = T)
groups = as(groups,'CsparseMatrix')
tree = list('eta_g'= eta_g,'groups' = groups,'own_variables' = own_variables,
  'N_own_variables' = N_own_variables)
.printf('\ntest prox tree-linf\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'tree-l2', pos = FALSE,intercept = FALSE)),n = 1)

.printf('\ntest prox tree-linf\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'tree-linf', pos = FALSE,intercept = FALSE)),n = 1)

.printf('Second tree example\n')
# Example 2 of tree structure
# tree structured groups:
# g1= {0 1 2 3 4 5 6 7 8 9}    root(g1) = { };
# g2= {0 1 2 3 4 5}            root(g2) = {0 1 2};
# g3= {3 4}                    root(g3) = {3 4};
# g4= {5}                      root(g4) = {5};
# g5= {6 7 8 9}                root(g5) = { };
# g6= {6 7}                    root(g6) = {6 7};
# g7= {8 9}                    root(g7) = {8};
# g8 = {9}                     root(g8) = {9};

own_variables = as.vector(c(0, 0, 3, 5, 6, 6, 8, 9),mode= 'integer')
N_own_variables = as.vector(c(0,3,2,1,0,2,1,1),mode= 'integer')
eta_g = as.vector(c(1,1,1,2,2,2,2.5,2.5),mode = 'double')
groups = matrix(as.vector(c(0,0,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,
  0,1,0,0,0,0,0,0,
  0,1,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,
  0,0,0,0,1,0,0,0,
  0,0,0,0,1,0,0,0,
  0,0,0,0,0,0,1,0),mode='logical'),ncol = 8,byrow = T)
groups = as(groups,'CsparseMatrix')
tree = list('eta_g'= eta_g,'groups' = groups,'own_variables' = own_variables,
  'N_own_variables' = N_own_variables)
.printf('\ntest prox tree-l0\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'tree-l0', pos = FALSE,intercept = FALSE)),n = 1)

.printf('\ntest prox tree-l2\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'tree-l2', pos = FALSE,intercept = FALSE)),n = 1)

.printf('\ntest prox tree-linf\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'tree-linf', pos = FALSE,intercept = FALSE)),n = 1)

# mexProximalTree also works with non-tree-structured regularization functions
.printf('\nprox l1, intercept, positivity constraint\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,regul = 'l1', pos = TRUE,intercept = TRUE)),n = 1)

.printf('\nprox multi-task tree\n')
alpha = Xtest1('spams',quote(spams.proximalTree(U,tree,FALSE,numThreads = -1,verbose = TRUE,lambda1 = 0.1,lambda2 = 0.1,regul = 'multi-task-tree', pos = FALSE,intercept = FALSE)),n = 1)

5.4  Function spams.proximalGraph

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

 
min
v ∈ ℝp
 
1
2
||uv||22 + λ 
 
g ∈ G
 ηg ||vg||,     (44)

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

 
min
V ∈ ℝp× n
 
1
2
||UV||F2 + λ 
n
i=1
 
g ∈ G
 ηg ||vgi||+ λ2
 
g ∈ G
 
 
max
j ∈ g
 ||vgj||,     (45)

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

#
# Name: spams.proximalGraph
#
# Usage: spams.proximalGraph(U,graph,return_val_loss = FALSE,numThreads =-1,lambda1=1.0,lambda2=0.,
#                     lambda3=0.,intercept=FALSE,resetflow=FALSE,regul="",verbose=FALSE,
#                     pos=FALSE,clever=TRUE,size_group=1,transpose=FALSE)
#
# Description:
#     spams.proximalGraph computes a proximal operator. Depending
#         on the value of regul, it computes 
#         
#         Given an input matrix U=[u^1,\ldots,u^n], and a set of groups G,
#         it computes a matrix V=[v^1,\ldots,v^n] such that
#         
#         if regul='graph'
#         for every column u of U, it computes a column v of V solving
#             argmin 0.5||u-v||_2^2 + lambda1\sum_{g \in G} \eta_g||v_g||_inf
#             
#         if regul='graph+ridge'
#         for every column u of U, it computes a column v of V solving
#             argmin 0.5||u-v||_2^2 + lambda1\sum_{g \in G} \eta_g||v_g||_inf + lambda1_2||v||_2^2
#             
#             
#         if regul='multi-task-graph'
#            V=argmin 0.5||U-V||_F^2 + lambda1 \sum_{i=1}^n\sum_{g \in G} \eta_g||v^i_g||_inf + ...
#                                                lambda1_2 \sum_{g \in G} \eta_g max_{j in g}||V_j||_{inf}
#                                                
#         it can also be used with any regularization addressed by spams.proximalFlat
#         
#         for all these regularizations, it is possible to enforce non-negativity constraints
#         with the option pos, and to prevent the last row of U to be regularized, with
#         the option intercept
#
# Inputs:
#       U:  double p x n matrix   (input signals)
#               m is the signal size
#       graph: struct
#             with three fields, eta_g, groups, and groups_var
#             
#             The first fields sets the weights for every group
#                graph.eta_g            double N vector 
#                
#             The next field sets inclusion relations between groups 
#             (but not between groups and variables):
#                graph.groups           sparse (double or boolean) N x N matrix  
#                the (i,j) entry is non-zero if and only if i is different than j and 
#                gi is included in gj.
#                
#             The next field sets inclusion relations between groups and variables
#                graph.groups_var       sparse (double or boolean) p x N matrix
#                the (i,j) entry is non-zero if and only if the variable i is included 
#                in gj, but not in any children of gj.
#                
#             examples are given in test_ProximalGraph.m
#             
#       return_val_loss:     
#               if true the function will return 2 matrices in a list.
#       lambda1:  (regularization parameter)
#       regul: (choice of regularization, see above)
#       lambda2:  (optional, regularization parameter)
#       lambda3:  (optional, regularization parameter)
#       verbose: (optional, verbosity level, false by default)
#       intercept: (optional, last row of U is not regularized,
#         false by default)
#       pos: (optional, adds positivity constraints on the
#         coefficients, false by default)
#       numThreads: (optional, number of threads for exploiting
#         multi-core / multi-cpus. By default, it takes the value -1,
#         which automatically selects all the available CPUs/cores).
#       resetflow:    undocumented; modify at your own risks!
#       clever:    undocumented; modify at your own risks!
#       size_group:    undocumented; modify at your own risks!
#       transpose:    undocumented; modify at your own risks!
#
# Output:
#       V: double p x n matrix (output coefficients)
#       val_regularizer: double 1 x n vector (value of the regularization
#       term at the optimum).
#       val_loss:        vector of size ncol(U)
#             alpha <- spams.proximalGraph(U,graph,return_val_loss = FALSE, ...)
#             v <- spams.proximalGraph(U,graph,return_val_loss = TRUE, ...)
#             alpha <- v[[1]]
#             val_loss <- v[[2]]
#
# Authors:
# Julien MAIRAL, 2010 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (R interface)
#
# Note:
#     Valid values for the regularization parameter (regul) are:
#       "l0", "l1", "l2", "linf", "l2-not-squared", "elastic-net", "fused-lasso",
#       "group-lasso-l2", "group-lasso-linf", "sparse-group-lasso-l2",
#       "sparse-group-lasso-linf", "l1l2", "l1linf", "l1l2+l1", "l1linf+l1",
#       "tree-l0", "tree-l2", "tree-linf", "graph", "graph-ridge", "graph-l2",
#       "multi-task-tree", "multi-task-graph", "l1linf-row-column", "trace-norm",
#       "trace-norm-vec", "rank", "rank-vec", "none"
#

The following piece of code illustrates how to use this function. The following piece of code contains usage examples:

library(spams)
set.seed(0)
lambda1 = 0.1 # regularization parameter
num_threads = -1 # all cores (-1 by default)
verbose = TRUE   # verbosity, false by default
pos = FALSE       # can be used with all the other regularizations
intercept = FALSE # can be used with all the other regularizations     

U = matrix(rnorm(1000),nrow = 10,ncol = 100,byrow = FALSE)
.printf('First graph example\n')
# Example 1 of graph structure
# groups:
# g1= {0 1 2 3}
# g2= {3 4 5 6}
# g3= {6 7 8 9}

eta_g = as.vector(c(1, 1, 1),mode='double')
groups = as(matrix(as.vector(rep(0,9),mode='logical'),nrow = 3,ncol = 3),'CsparseMatrix')
groups_var = as(matrix(as.vector(c(1, 0, 0,
  1, 0, 0,
  1, 0, 0,
  1, 1, 0,
  0, 1, 0,
  0, 1, 0,
  0, 1, 1,
  0, 0, 1,
  0, 0, 1,
  0, 0, 1),mode='logical'),ncol = 3,byrow = T),'CsparseMatrix')
graph = list('eta_g'= eta_g,'groups' = groups,'groups_var' = groups_var)
.printf('\ntest prox graph\n')
regul='graph'
alpha = Xtest1('spams',quote(spams.proximalGraph(U,graph,FALSE,lambda1 = lambda1,numThreads  = num_threads ,verbose = verbose,pos = pos,intercept = intercept,regul = regul)),n = 1)

# Example 2 of graph structure
# groups:
# g1= {0 1 2 3}
# g2= {3 4 5 6}
# g3= {6 7 8 9}
# g4= {0 1 2 3 4 5}
# g5= {6 7 8}


eta_g = as.vector(c(1, 1, 1, 1, 1),mode='double')
groups = as(matrix(as.vector(c(0, 0, 0, 1, 0,
  0, 0, 0, 0, 0,
  0, 0, 0, 0, 0,
  0, 0, 0, 0, 0,
  0, 0, 1, 0, 0),mode='logical'),ncol = 5,byrow = T),'CsparseMatrix')

groups_var = as(matrix(as.vector(c(1, 0, 0, 0, 0,
  1, 0, 0, 0, 0,
  1, 0, 0, 0, 0,
  1, 1, 0, 0, 0,
  0, 1, 0, 1, 0,
  0, 1, 0, 1, 0,
  0, 1, 0, 0, 1,
  0, 0, 0, 0, 1,
  0, 0, 0, 0, 1,
  0, 0, 1, 0, 0),mode='logical'),ncol = 5,byrow = T),'CsparseMatrix')
graph = list('eta_g'= eta_g,'groups' = groups,'groups_var' = groups_var)

.printf('\ntest prox graph\n')
alpha = Xtest1('spams',quote(spams.proximalGraph(U,graph,FALSE,lambda1 = lambda1,numThreads  = num_threads ,verbose = verbose,pos = pos,intercept = intercept,regul = regul)),n = 1)
 #
.printf('\ntest prox multi-task-graph\n')
regul = 'multi-task-graph'
lambda2 = 0.1
alpha = Xtest1('spams',quote(spams.proximalGraph(U,graph,FALSE,lambda1 = lambda1,lambda2 = lambda2,numThreads  = num_threads ,verbose = verbose,pos = pos,intercept = intercept,regul = regul)),n = 1)
 #
.printf('\ntest no regularization\n')
regul = 'none'
alpha = Xtest1('spams',quote(spams.proximalGraph(U,graph,FALSE,lambda1 = lambda1,lambda2 = lambda2,numThreads  = num_threads ,verbose = verbose,pos = pos,intercept = intercept,regul = regul)),n = 1)

5.5  Function spams.proximalPathCoding

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

#
# The R function is not yet implemented.
#

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 spams.evalPathCoding

#
# The R function is not yet implemented.
#

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 {xi,yi}i=1n, with xi ∈ ℝp and yi ∈ ℝ for all i in [ 1;n ], we address

 
min
w ∈ ℝpb ∈ ℝ
 
n
i=1
 
1
2
(yiwxib)2 + λψ(w),

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 ℝn that carries the entries yi. The problem without the intercept takes the following form, which we have already encountered in the previous toolbox, but with different notations:

 
min
w ∈ ℝp
 
1
2
||yXw||22 + λψ(w),

where the X=[xi,…,xn]T (the xi’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 {xi,yi}i=1n, with xi ∈ ℝp, but the variables yi are now in {−1,+1} for all i in [ 1;n ]. The optimization problem we address is

   
 
min
w ∈ ℝpb ∈ ℝ
 
1
n
n
i=1
 log(1+eyi(wxi+b) + λψ(w),

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 {xi,yi}i=1n, where the variables xi are still vectors in ℝp, but the yi’s have integer values in {1,2,…,r}. The formulation we address is the following multi-class learning problem

 
min
W ∈ ℝp × rb ∈ ℝr
 
1
n
n
i=1
 log

r
j=1
 e (wjwyi)xi + bjbyi

+ λ
r
j=1
ψ(wj),     (46)

where W = [w1,…,wr] 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 {xi,yi}i=1n, where the variables xi are still vectors in ℝp, and yi is a vector in ℝr. We are looking for r regression vectors wj, for j∈ [ 1;r ], or equivalently for a matrix W=[w1,…,wr] in ℝp × r. The formulation we address is the following multi-task regression problem

         
 
min
W ∈ ℝp × rb ∈ ℝr
 
r
j=1
n
i=1
 
1
2
(yjiwxibj)2 + λψ(W),

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

   
 
min
W ∈ ℝp × r
 
1
2
 ||YXW||F2 + λψ(W).

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 {xi,yi}i=1n, with the xi’s in ℝp, and the yi’s are vectors in {−1,+1}r. We look for a matrix W=[w1,…,wr] in ℝp × r. The formulation is the following multi-task regression problem

   
 
min
W ∈ ℝp × rb ∈ ℝr
 
r
j=1
1
n
n
i=1
 log

1+eyji(wxi+bj)

+ λψ(W).

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:

   
 
min
W ∈ ℝp × rb ∈ ℝr
 
1
n
n
i=1
 log

r
j=1
 e (wjwyi)xi + bjbyi

+ λψ(W).

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 spams.fistaFlat

Given a matrix X=[x1,…,xp]T in ℝm × p, and a matrix Y=[y1,…,yn], it solves the optimization problems presented in the previous section, with the same regularization functions as spams.proximalFlat. see usage details below:

#
# Name: spams.fistaFlat
#
# Usage: spams.fistaFlat(Y,X,W0,return_optim_info = FALSE,numThreads =-1,max_it =1000,L0=1.0,
#                 fixed_step=FALSE,gamma=1.5,lambda1=1.0,delta=1.0,lambda2=0.,lambda3=0.,
#                 a=1.0,b=0.,c=1.0,tol=0.000001,it0=100,max_iter_backtracking=1000,
#                 compute_gram=FALSE,lin_admm=FALSE,admm=FALSE,intercept=FALSE,
#                 resetflow=FALSE,regul="",loss="",verbose=FALSE,pos=FALSE,clever=FALSE,
#                 log=FALSE,ista=FALSE,subgrad=FALSE,logName="",is_inner_weights=FALSE,
#                 inner_weights=c(0.),size_group=1,groups = NULL,sqrt_step=TRUE,
#                 transpose=FALSE,linesearch_mode=0)
#
# Description:
#     spams.fistaFlat solves sparse regularized problems.
#         X is a design matrix of size m x p
#         X=[x^1,...,x^n]', where the x_i's are the rows of X
#         Y=[y^1,...,y^n] is a matrix of size m x n
#         It implements the algorithms FISTA, ISTA and subgradient descent.
#         
#           - if loss='square' and regul is a regularization function for vectors,
#             the entries of Y are real-valued,  W = [w^1,...,w^n] is a matrix of size p x n
#             For all column y of Y, it computes a column w of W such that
#               w = argmin 0.5||y- X w||_2^2 + lambda1 psi(w)
#               
#           - if loss='square' and regul is a regularization function for matrices
#             the entries of Y are real-valued,  W is a matrix of size p x n. 
#             It computes the matrix W such that
#               W = argmin 0.5||Y- X W||_F^2 + lambda1 psi(W)
#               
#           - loss='square-missing' same as loss='square', but handles missing data
#             represented by NaN (not a number) in the matrix Y
#             
#           - if loss='logistic' and regul is a regularization function for vectors,
#             the entries of Y are either -1 or +1, W = [w^1,...,w^n] is a matrix of size p x n
#             For all column y of Y, it computes a column w of W such that
#               w = argmin (1/m)sum_{j=1}^m log(1+e^(-y_j x^j' w)) + lambda1 psi(w),
#             where x^j is the j-th row of X.
#             
#           - if loss='logistic' and regul is a regularization function for matrices
#             the entries of Y are either -1 or +1, W is a matrix of size p x n
#               W = argmin sum_{i=1}^n(1/m)sum_{j=1}^m log(1+e^(-y^i_j x^j' w^i)) + lambda1 psi(W)
#               
#           - if loss='multi-logistic' and regul is a regularization function for vectors,
#             the entries of Y are in {0,1,...,N} where N is the total number of classes
#             W = [W^1,...,W^n] is a matrix of size p x Nn, each submatrix W^i is of size p x N
#             for all submatrix WW of W, and column y of Y, it computes
#               WW = argmin (1/m)sum_{j=1}^m log(sum_{j=1}^r e^(x^j'(ww^j-ww^{y_j}))) + lambda1 sum_{j=1}^N psi(ww^j),
#             where ww^j is the j-th column of WW.
#             
#           - if loss='multi-logistic' and regul is a regularization function for matrices,
#             the entries of Y are in {0,1,...,N} where N is the total number of classes
#             W is a matrix of size p x N, it computes
#               W = argmin (1/m)sum_{j=1}^m log(sum_{j=1}^r e^(x^j'(w^j-w^{y_j}))) + lambda1 psi(W)
#             where ww^j is the j-th column of WW.
#             
#           - loss='cur' useful to perform sparse CUR matrix decompositions, 
#               W = argmin 0.5||Y-X*W*X||_F^2 + lambda1 psi(W)
#               
#               
#         The function psi are those used by spams.proximalFlat (see documentation)
#         
#         This function can also handle intercepts (last row of W is not regularized),
#         and/or non-negativity constraints on W, and sparse matrices for X
#
# Inputs:
#       Y:  double dense m x n matrix
#       X:  double dense or sparse m x p matrix   
#       W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
#            initial guess
#       return_optim_info:     
#               if true the function will return 2 matrices in a list.
#       loss: (choice of loss, see above)
#       regul: (choice of regularization, see function spams.proximalFlat)
#       lambda1: (regularization parameter)
#       lambda2: (optional, regularization parameter, 0 by default)
#       lambda3: (optional, regularization parameter, 0 by default)
#       verbose: (optional, verbosity level, false by default)
#       pos: (optional, adds positivity constraints on the
#           coefficients, false by default)
#       transpose: (optional, transpose the matrix in the regularization function)
#       size_group: (optional, for regularization functions assuming a group
#            structure)
#       groups: (int32, optional, for regularization functions assuming a group
#            structure, see spams.proximalFlat)
#       numThreads: (optional, number of threads for exploiting
#           multi-core / multi-cpus. By default, it takes the value -1,
#           which automatically selects all the available CPUs/cores).
#       max_it: (optional, maximum number of iterations, 100 by default)
#       it0: (optional, frequency for computing duality gap, every 10 iterations by default)
#       tol: (optional, tolerance for stopping criteration, which is a relative duality gap
#           if it is available, or a relative change of parameters).
#       gamma: (optional, multiplier for increasing the parameter L in fista, 1.5 by default)
#       L0: (optional, initial parameter L in fista, 0.1 by default, should be small enough)
#       fixed_step: (deactive the line search for L in fista and use L0 instead)
#       linesearch_mode: (line-search scheme when ista=true:
#       0: default, monotonic backtracking scheme
#       1: monotonic backtracking scheme, with restart at each iteration
#       2: Barzilai-Borwein step sizes (similar to SparSA by Wright et al.) 
#       3: non-monotonic backtracking
#       compute_gram: (optional, pre-compute X^TX, false by default).
#       intercept: (optional, do not regularize last row of W, false by default).
#       ista: (optional, use ista instead of fista, false by default).
#       subgrad: (optional, if not ista, use subradient descent instead of fista, false by default).
#       a: 
#       b: (optional, if subgrad, the gradient step is a/(t+b)
#       also similar options as spams.proximalFlat
#       
#       the function also implements the ADMM algorithm via an option admm=true. It is not documented
#       and you need to look at the source code to use it.
#       delta:    undocumented; modify at your own risks!
#       c:    undocumented; modify at your own risks!
#       max_iter_backtracking:    undocumented; modify at your own risks!
#       lin_admm:    undocumented; modify at your own risks!
#       admm:    undocumented; modify at your own risks!
#       resetflow:    undocumented; modify at your own risks!
#       clever:    undocumented; modify at your own risks!
#       log:    undocumented; modify at your own risks!
#       logName:    undocumented; modify at your own risks!
#       is_inner_weights:    undocumented; modify at your own risks!
#       inner_weights:    undocumented; modify at your own risks!
#       sqrt_step:    undocumented; modify at your own risks!
#
# Output:
#       W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
#       optim: optional, double dense 4 x n matrix.
#           first row: values of the objective functions.
#           third row: values of the relative duality gap (if available)
#           fourth row: number of iterations
#       optim_info:        vector of size 4, containing information of the optimization.
#             W <- spams.fistaFlat(Y,X,W0,return_optim_info = FALSE,...)
#             v <- spams.fistaFlat(Y,X,W0,return_optim_info = TRUE,...)
#             W <- v[[1]]
#             optim_info <- v[[2]]
#
# Authors:
# Julien MAIRAL, 2010 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (R interface)
#
# Note:
#     Valid values for the regularization parameter (regul) are:
#       "l0", "l1", "l2", "linf", "l2-not-squared", "elastic-net", "fused-lasso",
#       "group-lasso-l2", "group-lasso-linf", "sparse-group-lasso-l2",
#       "sparse-group-lasso-linf", "l1l2", "l1linf", "l1l2+l1", "l1linf+l1",
#       "tree-l0", "tree-l2", "tree-linf", "graph", "graph-ridge", "graph-l2",
#       "multi-task-tree", "multi-task-graph", "l1linf-row-column", "trace-norm",
#       "trace-norm-vec", "rank", "rank-vec", "none"
#

The following piece of code illustrates how to use this function. The following piece of code contains usage examples:

library(spams)
set.seed(0)
m = 100;n = 200
X = matrix(rnorm(m * n),nrow = m,ncol = n,byrow = FALSE)
X = X - matrix(rep(colMeans(X),nrow(X)),nrow(X),ncol(X),byrow = T)
X = spams.normalize(X)
Y = matrix(rnorm(m),nrow = m,ncol = 1,byrow = FALSE)
Y = Y - matrix(rep(colMeans(Y),nrow(Y)),nrow(Y),ncol(Y),byrow = T)
Y = spams.normalize(Y)
W0 = matrix(c(0),nrow = ncol(X), ncol = ncol(Y))
# Regression experiments 
# 100 regression problems with the same design matrix X.

.printf("\nVarious regression experiments\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l1')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
 ###

.printf("\nFISTA + Regression l1\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l1',ista = TRUE)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
###
.printf("\nSubgradient Descent + Regression l1\n")
max_it = 200
it0 = 10
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = 50, max_it = 500,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l1',ista = FALSE,subgrad = TRUE,a = 0.1, b = 1000)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
###
.printf("\nFISTA + Regression l2\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = 50, max_it = 500,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l2',ista = FALSE,subgrad = TRUE,a = 0.1, b = 1000)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
###
.printf("\nFISTA + Regression l2 + sparse feature matrix\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,as(X,'CsparseMatrix'),W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l2',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
#######
.printf("\nFISTA + Regression Elastic-Net\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'elastic-net',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

.printf("\nFISTA + Group Lasso L2\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'group-lasso-l2',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,size_group = 2)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

.printf("\nFISTA + Group Lasso L2 with variable size of groups\n")
groups = matrix(sample(1:5,ncol(X),replace=TRUE),nrow = 1)
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.5, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'group-lasso-l2',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,size_group = 2,
  groups = groups)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])
.printf("\nFISTA + Trace Norm\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'trace-norm-vec',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

###
.printf("\nFISTA + Regression Fused-Lasso\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'fused-lasso',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

.printf("\nFISTA + Regression no regularization\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'none',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

.printf("\nFISTA + Regression l1 with intercept\n")
x1 = cbind(X,matrix(1,nrow = nrow(X),ncol = 1))
W01 = rbind(W0,matrix(0,nrow = 1, ncol = ncol(W0)))
res = Xtest1('spams',quote(spams.fistaFlat(Y,x1,W01,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = TRUE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

.printf("\nFISTA + Regression l1 with intercept+ non-negative\n")
x1 = cbind(X,matrix(1,nrow = nrow(X),ncol = 1))
W01 = rbind(W0,matrix(0,nrow = 1, ncol = ncol(W0)))
res = Xtest1('spams',quote(spams.fistaFlat(Y,x1,W01,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = TRUE,pos = TRUE,compute_gram = TRUE, loss = 'square',regul = 'l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

.printf("\nISTA + Regression l0\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.05, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'l0',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])
# Classification

.printf("\nOne classification experiment\n")
Y = matrix(2. * as.double(rnorm(100) > 0.) - 1.,nrow = 100,ncol = 1,byrow = FALSE)
.printf("\nFISTA + Logistic l1\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'logistic',regul = 'l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# can be used of course with other regularization functions, intercept,...
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'weighted-logistic',regul = 'l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# can be used of course with other regularization functions, intercept,...

.printf("\nFISTA + Logistic l1 + sparse matrix\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,as(X,'CsparseMatrix'),W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'logistic',regul = 'l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# can be used of course with other regularization functions, intercept,...

# Multi-Class classification

Y = ceiling(5 * matrix(runif(100 * 1000,0,1),nrow = 100,ncol = 1000,byrow = FALSE)) - 1
.printf("\nFISTA + Multi-Class Logistic l1\n")
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'multi-logistic',regul = 'l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# can be used of course with other regularization functions, intercept,...

# Multi-Task regression

Y = matrix(rnorm(100 * 100),nrow = 100,ncol = 100,byrow = FALSE)
Y = Y - matrix(rep(colMeans(Y),nrow(Y)),nrow(Y),ncol(Y),byrow = T)
Y = spams.normalize(Y)
W0 = matrix(c(0),nrow = ncol(X), ncol = ncol(Y))
.printf("\nFISTA + Regression l1l2\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'square',regul = 'l1l2',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

.printf("\nFISTA + Regression l1linf\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'square',regul = 'l1linf',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])


.printf("\nFISTA + Regression l1l2 + l1\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'square',regul = 'l1l2+l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

.printf("\nFISTA + Regression l1linf + l1\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'square',regul = 'l1linf+l1',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

.printf("\nFISTA + Regression l1linf + row + columns\n")
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'square',regul = 'l1linf-row-column',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

 # Multi-Task Classification

.printf("\nFISTA + Logistic + l1l2\n")
Y = matrix(2. * as.double(rnorm(100 * 100) > 1.) - 1.,nrow = 100,ncol = 100,byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'logistic',regul = 'l1l2',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# Multi-Class + Multi-Task Regularization

.printf("\nFISTA + Multi-Class Logistic l1l2\n")
Y = ceiling(5 * matrix(runif(100 * 1000,0,1),nrow = 100,ncol = 1000,byrow = FALSE)) - 1
Y = spams.normalize(Y)
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaFlat(Y,X,W0,TRUE,numThreads = 1,verbose = TRUE,lambda1 = 0.01, it0 = it0, max_it = max_it,L0 = 0.1, tol = 1e-3, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'multi-logistic',regul = 'l1l2',ista = FALSE,subgrad = FALSE,a = 0.1, b = 1000,lambda2 = 0.1,lambda3 = 0.1,size_group = 5)),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# can be used of course with other regularization functions, intercept,...

5.9  Function spams.fistaTree

Given a matrix X=[x1,…,xp]T in ℝm × p, and a matrix Y=[y1,…,yn], it solves the optimization problems presented in the previous section, with the same regularization functions as spams.proximalTree. see usage details below:

#
# Name: spams.fistaTree
#
# Usage: spams.fistaTree(Y,X,W0,tree,return_optim_info = FALSE,numThreads =-1,max_it =1000,L0=1.0,
#                 fixed_step=FALSE,gamma=1.5,lambda1=1.0,delta=1.0,lambda2=0.,lambda3=0.,
#                 a=1.0,b=0.,c=1.0,tol=0.000001,it0=100,max_iter_backtracking=1000,
#                 compute_gram=FALSE,lin_admm=FALSE,admm=FALSE,intercept=FALSE,
#                 resetflow=FALSE,regul="",loss="",verbose=FALSE,pos=FALSE,clever=FALSE,
#                 log=FALSE,ista=FALSE,subgrad=FALSE,logName="",is_inner_weights=FALSE,
#                 inner_weights=c(0.),size_group=1,sqrt_step=TRUE,transpose=FALSE,
#                 linesearch_mode=0)
#
# Description:
#     spams.fistaTree solves sparse regularized problems.
#         X is a design matrix of size m x p
#         X=[x^1,...,x^n]', where the x_i's are the rows of X
#         Y=[y^1,...,y^n] is a matrix of size m x n
#         It implements the algorithms FISTA, ISTA and subgradient descent for solving
#         
#           min_W  loss(W) + lambda1 psi(W)
#           
#         The function psi are those used by spams.proximalTree (see documentation)
#         for the loss functions, see the documentation of spams.fistaFlat
#         
#         This function can also handle intercepts (last row of W is not regularized),
#         and/or non-negativity constraints on W and sparse matrices X
#
# Inputs:
#       Y:  double dense m x n matrix
#       X:  double dense or sparse m x p matrix   
#       W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
#            initial guess
#       tree: named list (see documentation of spams.proximalTree)
#       return_optim_info:     
#               if true the function will return 2 matrices in a list.
#       loss: (choice of loss, see above)
#       regul: (choice of regularization, see function spams.proximalFlat)
#       lambda1: (regularization parameter)
#       lambda2: (optional, regularization parameter, 0 by default)
#       lambda3: (optional, regularization parameter, 0 by default)
#       verbose: (optional, verbosity level, false by default)
#       pos: (optional, adds positivity constraints on the
#           coefficients, false by default)
#       transpose: (optional, transpose the matrix in the regularization function)
#       size_group: (optional, for regularization functions assuming a group
#            structure)
#       numThreads: (optional, number of threads for exploiting
#           multi-core / multi-cpus. By default, it takes the value -1,
#           which automatically selects all the available CPUs/cores).
#       max_it: (optional, maximum number of iterations, 100 by default)
#       it0: (optional, frequency for computing duality gap, every 10 iterations by default)
#       tol: (optional, tolerance for stopping criteration, which is a relative duality gap
#           if it is available, or a relative change of parameters).
#       gamma: (optional, multiplier for increasing the parameter L in fista, 1.5 by default)
#       L0: (optional, initial parameter L in fista, 0.1 by default, should be small enough)
#       fixed_step: (deactive the line search for L in fista and use L0 instead)
#       compute_gram: (optional, pre-compute X^TX, false by default).
#       intercept: (optional, do not regularize last row of W, false by default).
#       ista: (optional, use ista instead of fista, false by default).
#       subgrad: (optional, if not ista, use subradient descent instead of fista, false by default).
#       a: 
#       b: (optional, if subgrad, the gradient step is a/(t+b)
#       also similar options as spams.proximalTree
#       
#       the function also implements the ADMM algorithm via an option admm=true. It is not documented
#       and you need to look at the source code to use it.
#       delta:    undocumented; modify at your own risks!
#       c:    undocumented; modify at your own risks!
#       max_iter_backtracking:    undocumented; modify at your own risks!
#       lin_admm:    undocumented; modify at your own risks!
#       admm:    undocumented; modify at your own risks!
#       resetflow:    undocumented; modify at your own risks!
#       clever:    undocumented; modify at your own risks!
#       log:    undocumented; modify at your own risks!
#       logName:    undocumented; modify at your own risks!
#       is_inner_weights:    undocumented; modify at your own risks!
#       inner_weights:    undocumented; modify at your own risks!
#       sqrt_step:    undocumented; modify at your own risks!
#
# Output:
#       W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
#       optim: optional, double dense 4 x n matrix.
#           first row: values of the objective functions.
#           third row: values of the relative duality gap (if available)
#           fourth row: number of iterations
#       optim_info:        vector of size 4, containing information of the optimization.
#             W <- spams.fistaTree(Y,X,W0,tree,return_optim_info = FALSE,...)
#             v <- spams.fistaTree(Y,X,W0,tree,return_optim_info = TRUE,...)
#             W <- v[[1]]
#             optim_info <- v[[2]]
#
# Authors:
# Julien MAIRAL, 2010 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (R interface)
#
# Note:
#     Valid values for the regularization parameter (regul) are:
#       "l0", "l1", "l2", "linf", "l2-not-squared", "elastic-net", "fused-lasso",
#       "group-lasso-l2", "group-lasso-linf", "sparse-group-lasso-l2",
#       "sparse-group-lasso-linf", "l1l2", "l1linf", "l1l2+l1", "l1linf+l1",
#       "tree-l0", "tree-l2", "tree-linf", "graph", "graph-ridge", "graph-l2",
#       "multi-task-tree", "multi-task-graph", "l1linf-row-column", "trace-norm",
#       "trace-norm-vec", "rank", "rank-vec", "none"
#

The following piece of code illustrates how to use this function. The following piece of code contains usage examples:

library(spams)
set.seed(0)
m = 100;n = 10
X = matrix(rnorm(m * n),nrow = m,ncol = n,byrow = FALSE)
X = X - matrix(rep(colMeans(X),nrow(X)),nrow(X),ncol(X),byrow = T)
X = spams.normalize(X)
Y = matrix(rnorm(m),nrow = m,ncol = 1,byrow = FALSE)
Y = Y - matrix(rep(colMeans(Y),nrow(Y)),nrow(Y),ncol(Y),byrow = T)
Y = spams.normalize(Y)
W0 = matrix(c(0),nrow = ncol(X), ncol = ncol(Y))
own_variables = as.vector(c(0,0,3,5,6,6,8,9),mode= 'integer')
N_own_variables = as.vector(c(0,3,2,1,0,2,1,1),mode= 'integer')
eta_g = as.vector(c(1,1,1,2,2,2,2.5,2.5),mode = 'double')
groups = matrix(as.vector(c(0, 0, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 0, 0, 0, 0, 0,
  0, 1, 0, 0, 0, 0, 0, 0,
  0, 1, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 1, 0, 0, 0,
  0, 0, 0, 0, 1, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 1, 0),mode='logical'),ncol = 8,byrow = T)
groups = as(groups,'CsparseMatrix')
tree = list('eta_g'= eta_g,'groups' = groups,'own_variables' = own_variables,
          'N_own_variables' = N_own_variables)
.printf('\nVarious regression experiments\n')
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'tree-l2')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])
 ###
.printf('\nFISTA + Regression tree-linf\n')
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'tree-linf')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
###
# works also with non tree-structured regularization. tree is ignored

.printf('\nFISTA + Regression Fused-Lasso\n')
 res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'fused-lasso')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])
 ###
.printf('\nISTA + Regression tree-l0\n')
 res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'tree-l0')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])
 ###
.printf('\nFISTA + Regression tree-l2 with intercept\n')
x1 = cbind(X,matrix(1,nrow = nrow(X),ncol = 1))
W01 = rbind(W0,matrix(0,nrow = 1, ncol = ncol(W0)))
res = Xtest1('spams',quote(spams.fistaTree(Y,x1,W01,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = TRUE,pos = FALSE,compute_gram = TRUE, loss = 'square',regul = 'tree-l2')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])

#    Classification

.printf('\nOne classification experiment')
Y = matrix(2. * as.double(rnorm(100 * ncol(Y)) > 0.) - 1.,nrow = 100,ncol = ncol(Y),byrow = FALSE)
.printf('\nFISTA + Logistic + tree-linf\n')
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'logistic',regul = 'tree-linf')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
###
# can be used of course with other regularization functions, intercept,...

#  Multi-Class classification

Y = ceiling(5 * matrix(runif(100 * ncol(Y),0,1),nrow = 100,ncol = ncol(Y),byrow = FALSE)) - 1
.printf('\nFISTA + Multi-Class Logistic + tree-l2\n')
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = TRUE, loss = 'multi-logistic',regul = 'tree-l2')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, number of iterations: %f\n",optim_info[1],optim_info[4])
# can be used of course with other regularization functions, intercept,...

# Multi-Task regression

Y = matrix(rnorm(100 * 100),nrow = 100,ncol = 100,byrow = FALSE)
Y = Y - matrix(rep(colMeans(Y),nrow(Y)),nrow(Y),ncol(Y),byrow = T)
Y = spams.normalize(Y)
W0 = matrix(c(0),nrow = ncol(X), ncol = ncol(Y))
.printf('\nFISTA + Regression  multi-task-tree\n')
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = TRUE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'square',regul = 'multi-task-tree')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

# Multi-Task Classification
.printf('\nFISTA + Logistic + multi-task-tree\n')
Y = matrix(rnorm(100 * ncol(Y)),nrow = 100,ncol = ncol(Y),byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = TRUE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'logistic',regul = 'multi-task-tree')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

#  Multi-Class + Multi-Task Regularization
.printf('\nFISTA + Multi-Class Logistic +multi-task-tree\n')
Y = ceiling(5 * matrix(runif(100 * ncol(Y),0,1),nrow = 100,ncol = ncol(Y),byrow = FALSE)) - 1
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaTree(Y,X,W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'multi-logistic',regul = 'multi-task-tree')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])
# can be used of course with other regularization functions, intercept,...

.printf('\nFISTA + Multi-Class Logistic +multi-task-tree + sparse matrix\n')
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
res = Xtest1('spams',quote(spams.fistaTree(Y,as(X,'CsparseMatrix'),W0,tree,TRUE,numThreads = -1,verbose = FALSE,lambda1 = 0.001,lambda2 = 0.001,lambda3 = 0.001, it0 = 10, max_it = 200,L0 = 0.1, tol = 1e-5, intercept = FALSE,pos = FALSE,compute_gram = FALSE, loss = 'multi-logistic',regul = 'multi-task-tree')),n = 1)
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, number of iterations: %f\n",optim_info[1],optim_info[3],optim_info[4])

5.10  Function spams.fistaGraph

Given a matrix X=[x1,…,xp]T in ℝm × p, and a matrix Y=[y1,…,yn], it solves the optimization problems presented in the previous section, with the same regularization functions as spams.proximalGraph. see usage details below:

#
# Name: spams.fistaGraph
#
# Usage: spams.fistaGraph(Y,X,W0,graph,return_optim_info = FALSE,numThreads =-1,max_it =1000,L0=1.0,
#                  fixed_step=FALSE,gamma=1.5,lambda1=1.0,delta=1.0,lambda2=0.,lambda3=0.,
#                  a=1.0,b=0.,c=1.0,tol=0.000001,it0=100,max_iter_backtracking=1000,
#                  compute_gram=FALSE,lin_admm=FALSE,admm=FALSE,intercept=FALSE,
#                  resetflow=FALSE,regul="",loss="",verbose=FALSE,pos=FALSE,clever=FALSE,
#                  log=FALSE,ista=FALSE,subgrad=FALSE,logName="",is_inner_weights=FALSE,
#                  inner_weights=c(0.),size_group=1,sqrt_step=TRUE,transpose=FALSE,
#                  linesearch_mode=0)
#
# Description:
#     spams.fistaGraph solves sparse regularized problems.
#         X is a design matrix of size m x p
#         X=[x^1,...,x^n]', where the x_i's are the rows of X
#         Y=[y^1,...,y^n] is a matrix of size m x n
#         It implements the algorithms FISTA, ISTA and subgradient descent.
#         
#         It implements the algorithms FISTA, ISTA and subgradient descent for solving
#         
#           min_W  loss(W) + lambda1 psi(W)
#           
#         The function psi are those used by spams.proximalGraph (see documentation)
#         for the loss functions, see the documentation of spams.fistaFlat
#         
#         This function can also handle intercepts (last row of W is not regularized),
#         and/or non-negativity constraints on W.
#
# Inputs:
#       Y:  double dense m x n matrix
#       X:  double dense or sparse m x p matrix   
#       W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
#            initial guess
#       graph: struct (see documentation of spams.proximalGraph)
#       return_optim_info:     
#               if true the function will return 2 matrices in a list.
#       loss: (choice of loss, see above)
#       regul: (choice of regularization, see function spams.proximalFlat)
#       lambda1: (regularization parameter)
#       lambda2: (optional, regularization parameter, 0 by default)
#       lambda3: (optional, regularization parameter, 0 by default)
#       verbose: (optional, verbosity level, false by default)
#       pos: (optional, adds positivity constraints on the
#           coefficients, false by default)
#       numThreads: (optional, number of threads for exploiting
#           multi-core / multi-cpus. By default, it takes the value -1,
#           which automatically selects all the available CPUs/cores).
#       max_it: (optional, maximum number of iterations, 100 by default)
#       it0: (optional, frequency for computing duality gap, every 10 iterations by default)
#       tol: (optional, tolerance for stopping criteration, which is a relative duality gap
#           if it is available, or a relative change of parameters).
#       gamma: (optional, multiplier for increasing the parameter L in fista, 1.5 by default)
#       L0: (optional, initial parameter L in fista, 0.1 by default, should be small enough)
#       fixed_step: (deactive the line search for L in fista and use L0 instead)
#       compute_gram: (optional, pre-compute X^TX, false by default).
#       intercept: (optional, do not regularize last row of W, false by default).
#       ista: (optional, use ista instead of fista, false by default).
#       subgrad: (optional, if not ista, use subradient descent instead of fista, false by default).
#       a: 
#       b: (optional, if subgrad, the gradient step is a/(t+b)
#       also similar options as spams.proximalTree
#       
#       the function also implements the ADMM algorithm via an option admm=true. It is not documented
#       and you need to look at the source code to use it.
#       delta:    undocumented; modify at your own risks!
#       c:    undocumented; modify at your own risks!
#       max_iter_backtracking:    undocumented; modify at your own risks!
#       lin_admm:    undocumented; modify at your own risks!
#       admm:    undocumented; modify at your own risks!
#       resetflow:    undocumented; modify at your own risks!
#       clever:    undocumented; modify at your own risks!
#       log:    undocumented; modify at your own risks!
#       logName:    undocumented; modify at your own risks!
#       is_inner_weights:    undocumented; modify at your own risks!
#       inner_weights:    undocumented; modify at your own risks!
#       sqrt_step:    undocumented; modify at your own risks!
#       size_group:    undocumented; modify at your own risks!
#       transpose:    undocumented; modify at your own risks!
#
# Output:
#       W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)
#       optim: optional, double dense 4 x n matrix.
#           first row: values of the objective functions.
#           third row: values of the relative duality gap (if available)
#           fourth row: number of iterations
#       optim_info:        vector of size 4, containing information of the optimization.
#             W <- spams.fistaGraph(Y,X,W0,graph,return_optim_info = FALSE,...)
#             v <- spams.fistaGraph(Y,X,W0,graph,return_optim_info = TRUE,...)
#             W <- v[[1]]
#             optim_info <- v[[2]]
#
# Authors:
# Julien MAIRAL, 2010 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (R interface)
#
# Note:
#     Valid values for the regularization parameter (regul) are:
#       "l0", "l1", "l2", "linf", "l2-not-squared", "elastic-net", "fused-lasso",
#       "group-lasso-l2", "group-lasso-linf", "sparse-group-lasso-l2",
#       "sparse-group-lasso-linf", "l1l2", "l1linf", "l1l2+l1", "l1linf+l1",
#       "tree-l0", "tree-l2", "tree-linf", "graph", "graph-ridge", "graph-l2",
#       "multi-task-tree", "multi-task-graph", "l1linf-row-column", "trace-norm",
#       "trace-norm-vec", "rank", "rank-vec", "none"
#

The following piece of code illustrates how to use this function. The following piece of code contains usage examples:

library(spams)
set.seed(0)
num_threads = -1 # all cores (-1 by default)
verbose = FALSE   # verbosity, false by default
lambda1 = 0.1 # regularization ter
it0 = 1      # frequency for duality gap computations
max_it = 100 # maximum number of iterations
L0 = 0.1
tol = 1e-5
intercept = FALSE
pos = FALSE
eta_g = as.vector(c(1, 1, 1, 1, 1),mode='double')
groups = as(matrix(as.vector(c(0, 0, 0, 1, 0,
  0, 0, 0, 0, 0,
  0, 0, 0, 0, 0,
  0, 0, 0, 0, 0,
  0, 0, 1, 0, 0),mode='logical'),ncol = 5,byrow = T),'CsparseMatrix')
groups_var = as(matrix(as.vector(c(1, 0, 0, 0, 0,
  1, 0, 0, 0, 0,
  1, 0, 0, 0, 0,
  1, 1, 0, 0, 0,
  0, 1, 0, 1, 0,
  0, 1, 0, 1, 0,
  0, 1, 0, 0, 1,
  0, 0, 0, 0, 1,
  0, 0, 0, 0, 1,
  0, 0, 1, 0, 0),mode='logical'),ncol = 5,byrow = T),'CsparseMatrix')
graph = list('eta_g'= eta_g,'groups' = groups,'groups_var' = groups_var)
verbose = TRUE
X = matrix(rnorm(1000),nrow = 100,ncol = 10,byrow = FALSE)
X = X - matrix(rep(colMeans(X),nrow(X)),nrow(X),ncol(X),byrow = T)
X = spams.normalize(X)
Y = matrix(rnorm(100),nrow = 100,ncol = 1,byrow = FALSE)
Y = Y - matrix(rep(colMeans(Y),nrow(Y)),nrow(Y),ncol(Y),byrow = T)
Y = spams.normalize(Y)
W0 = matrix(c(0),nrow = ncol(X), ncol = ncol(Y))
 # Regression experiments 
# 100 regression problems with the same design matrix X.

.printf('\nVarious regression experiments\n')
compute_gram = TRUE
#
.printf('\nFISTA + Regression graph\n')
loss = 'square'
regul = 'graph'
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
.printf('\nADMM + Regression graph\n')
admm = TRUE
lin_admm = TRUE
c = 1
delta = 1
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,delta = delta)
tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
admm = FALSE
max_it = 5
it0 = 1
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,delta = delta)
tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
#  works also with non graph-structured regularization. graph is ignored

.printf('\nFISTA + Regression Fused-Lasso\n')
regul = 'fused-lasso'
lambda2 = 0.01
lambda3 = 0.01
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)
tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f,  time: %f, number of iterations: %f\n",optim_info[1],t,optim_info[4])
#
.printf("\nFISTA + Regression graph with intercept\n")
regul = 'graph'
intercept = TRUE
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
intercept = FALSE
# Classification
.printf('\nOne classification experiment\n')
Y = matrix(2. * as.double(rnorm(100 * ncol(Y)) > 0.) - 1.,nrow = 100,ncol = ncol(Y),byrow = FALSE)
.printf('\nFISTA + Logistic + graph-linf\n')
loss = 'logistic'
regul = 'graph'
lambda1 = 0.01
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
# can be used of course with other regularization functions, intercept,...

# Multi-Class classification

Y = ceiling(5 * matrix(runif(100 * ncol(Y),0,1),nrow = 100,ncol = ncol(Y),byrow = FALSE)) - 1
loss = 'multi-logistic'
regul = 'graph'
.printf('\nFISTA + Multi-Class Logistic + graph\n')
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
# can be used of course with other regularization functions, intercept,...

# Multi-Task regression

Y = matrix(rnorm(100 * ncol(Y)),nrow = 100,ncol = ncol(Y),byrow = FALSE)
Y = Y - matrix(rep(colMeans(Y),nrow(Y)),nrow(Y),ncol(Y),byrow = T)
Y = spams.normalize(Y)
W0 = matrix(c(0),nrow = ncol(X), ncol = ncol(Y))
compute_gram = FALSE
verbose = TRUE
loss = 'square'
.printf('\nFISTA + Regression  multi-task-graph\n')
regul = 'multi-task-graph'
lambda2 = 0.01
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
# Multi-Task Classification

.printf('\nFISTA + Logistic + multi-task-graph\n')
regul = 'multi-task-graph'
lambda2 = 0.01
loss = 'logistic'
Y = matrix(rnorm(100 * ncol(Y)),nrow = 100,ncol = ncol(Y),byrow = FALSE)
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
# Multi-Class + Multi-Task Regularization

verbose = FALSE
.printf('\nFISTA + Multi-Class Logistic +multi-task-graph\n')
Y = ceiling(5 * matrix(runif(100 * ncol(Y),0,1),nrow = 100,ncol = ncol(Y),byrow = FALSE)) - 1
nclasses = max(Y) + 1
W0 = matrix(0,nrow = ncol(X),nclasses * ncol(Y),byrow = FALSE)
loss = 'multi-logistic'
regul = 'multi-task-graph'
tic = proc.time()
res = spams.fistaGraph(
  Y,X,W0,graph,TRUE,numThreads = num_threads,verbose = verbose,
  lambda1 = lambda1,it0 = it0,max_it = max_it,L0 = L0,tol = tol,
  intercept = intercept,pos = pos,compute_gram = compute_gram,
  loss = loss,regul = regul,admm = admm,lin_admm = lin_admm,
  c = c,lambda2 = lambda2,lambda3 = lambda3,delta = delta)

tac = proc.time()
t = (tac - tic)[['elapsed']]
W = res[[1]]
optim_info = res[[2]]
.printf("mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n",optim_info[1],optim_info[3],t,optim_info[4])
#
# can be used of course with other regularization functions, intercept,...

5.11  Function spams.fistaPathCoding

Similarly, the toolbox handles the penalties of [24] with the following function

#
# The R function is not yet implemented.
#

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

5.12  Function solverPoisson

The following problem is addressed here

   
 
min
w ∈ ℝ+p
  
p
j=1
 xiw + δ − yi log(xiw +δ) + λ ψ(w).
#
# The R function is not yet implemented.
#

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

5.13  Function mexIncrementalProx

This implements the incremental solver MISO [25].

#
# The R function is not yet implemented.
#

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

5.14  Function mexStochasticProx

This implements the stochastic proximal gradient solver [26].

#
# The R function is not yet implemented.
#

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


Previous Up Next