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.
This is a fast implementation of the Orthogonal Matching Pursuit algorithm (or forward selection) [27, 35]. Given a matrix of signals X=[x1,…,xn] in ℝm × n and a dictionary D=[d1,…,dp] in ℝm × p, the algorithm computes a matrix A=[α1,…,αn] 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
| ||x−Dα||22 s.t. ||α||0 ≤ L, (2) |
or
| ||α||0 s.t. ||x−Dα||22 ≤ ε, (3) |
or
|
| ||x−Dα||22 + λ ||α||0. (4) |
For efficienty reasons, the method first computes the covariance matrix DTD, then for each signal, it computes DTx and performs the decomposition with a Cholesky-based algorithm (see [6] for instance).
Note that spams.omp can return the “greedy” regularization path if needed (see below):
| # # Name: spams.omp # # Usage: spams.omp(X,D,L = NULL,eps = NULL,lambda1 = NULL,return_reg_path = FALSE,numThreads = -1) # # Description: # spams.omp is an efficient implementation of the # Orthogonal Matching Pursuit algorithm. It is optimized # for solving a large number of small or medium-sized # decomposition problem (and not for a single large one). # It first computes the Gram matrix D'D and then perform # a Cholesky-based OMP of the input signals in parallel. # X=[x^1,...,x^n] is a matrix of signals, and it returns # a matrix A=[alpha^1,...,alpha^n] of coefficients. # # it addresses for all columns x of X, # min_{alpha} ||alpha||_0 s.t ||x-Dalpha||_2^2 <= eps # or # min_{alpha} ||x-Dalpha||_2^2 s.t. ||alpha||_0 <= L # or # min_{alpha} 0.5||x-Dalpha||_2^2 + lambda1||alpha||_0 # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # All the columns of D should have unit-norm ! # return_reg_path: # if true the function will return 2 matrices in a list. # L: (optional, maximum number of elements in each decomposition, # min(m,p) by default) # eps: (optional, threshold on the squared l2-norm of the residual, # 0 by default # lambda1: (optional, penalty parameter, 0 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). # # Output: # A: double sparse p x n matrix (output coefficients) # path (optional): double dense p x L matrix (regularization path of the first signal) # A <- spams.omp(X,D,L,eps,return_reg_path = FALSE,...) # v <- spams.omp(X,D,L,eps,return_reg_path = TRUE,...) # A <- v[[1]] # path <- v[[2]] # # Authors: # Julien MAIRAL, 2009 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # - Passing an int32 vector of length n to L provides # a different parameter L for each input signal x_i # - Passing a double vector of length n to eps and or lambda1 # provides a different parameter eps (or lambda1) for each input signal x_i # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) .printf("test omp\n") X = matrix(rnorm(64 * 100000),nrow = 64,ncol = 100000,byrow = FALSE) D = matrix(rnorm(64 * 200),nrow = 64,ncol = 200,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) L = 10 eps =0.1 numThreads = -1 tic = proc.time() alpha = spams.omp(X,D,L=L,eps=eps,return_reg_path = FALSE,numThreads = numThreads) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) ######################################## # Regularization path of a single signal ######################################## X = matrix(rnorm(64 * 1),nrow = 64,ncol = 1,byrow = FALSE) D = matrix(rnorm(64 * 10),nrow = 64,ncol = 10,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) L = 5 res = spams.omp(X,D,L = L,eps = eps,return_reg_path = TRUE,numThreads = numThreads) alpha = res[[1]] path = res[[2]] |
This is a variant of spams.omp with the possibility of handling a binary mask. Given a binary mask B=[β1,…,βn] in {0,1}m × n, it returns a matrix A=[α1,…,αn] such that for every column x of X, β of B, it computes a column α of A by addressing
| ||diag(β)(x−Dα)||22 s.t. ||α||0 ≤ L, (5) |
or
| ||α||0 s.t. ||diag(β)(x−Dα)||22 ≤ ε |
| , (6) |
or
|
| ||diag(β)(x−Dα)||22 +λ||α||0. (7) |
where diag(β) is a diagonal matrix with the entries of β on the diagonal.
| # # Name: spams.ompMask # # Usage: spams.ompMask(X,D,B,L = NULL,eps = NULL,lambda1 = NULL,return_reg_path = FALSE, # numThreads = -1) # # Description: # spams.ompMask is a variant of spams.omp that allow using # a binary mask B # # for all columns x of X, and columns beta of B, it computes a column # alpha of A by addressing # min_{alpha} ||alpha||_0 s.t ||diag(beta)*(x-Dalpha)||_2^2 # <= eps*||beta||_0/m # or # min_{alpha} ||diag(beta)*(x-Dalpha)||_2^2 s.t. ||alpha||_0 <= L # or # min_{alpha} 0.5||diag(beta)*(x-Dalpha)||_2^2 + lambda1||alpha||_0 # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # All the columns of D should have unit-norm ! # B: boolean m x n matrix (mask) # p is the number of elements in the dictionary # return_reg_path: # if true the function will return 2 matrices in a list. # L: (optional, maximum number of elements in each decomposition, # min(m,p) by default) # eps: (optional, threshold on the squared l2-norm of the residual, # 0 by default # lambda1: (optional, penalty parameter, 0 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). # # Output: # A: double sparse p x n matrix (output coefficients) # path (optional): double dense p x L matrix # (regularization path of the first signal) # A <- spams.ompMask(X,D,B,L,eps,return_reg_path = FALSE,...) # v <- spams.ompMask(X,D,B,L,eps,return_reg_path = TRUE,...) # A <- v[[1]] # path <- v[[2]] # # Authors: # Julien MAIRAL, 2010 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # - Passing an int32 vector of length n to L provides # a different parameter L for each input signal x_i # - Passing a double vector of length n to eps and or lambda1 # provides a different parameter eps (or lambda1) for each input signal x_i # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) .printf("test ompMask\n") X = matrix(rnorm(100 * 100),nrow = 100,ncol = 100,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) D = matrix(rnorm(100 * 20),nrow = 100,ncol = 20,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) mask = (X > 0) # generating a binary mask L = 20 eps =0.01 numThreads = -1 tic = proc.time() alpha = spams.ompMask(X,D,mask,L = L,eps = eps,return_reg_path = FALSE,numThreads = numThreads) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) |
This is a ridge regression solver using a conjugate gradient solver.
| # # The R function is not yet implemented. # |
The following piece of code illustrates how to use this function.
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=[x1,…,xn] in ℝm × n and a dictionary D in ℝm × p, depending on the input parameters, the algorithm returns a matrix of coefficients A=[α1,…,αn] in ℝp × n such that for every column x of X, the corresponding column α of A is the solution of
| ||x−Dα||22 s.t. ||α||1 ≤ λ, (8) |
or
| ||α||1 s.t. ||x−Dα||22 ≤ λ, (9) |
or
|
| ||x−Dα||22 + λ ||α||1 + |
| ||α||22. (10) |
For efficiency reasons, the method first compute the covariance matrix DTD, then for each signal, it computes DTx 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 spams.lasso can return the whole regularization path of the first signal x1 and can handle implicitely the matrix D if the quantities DTD and DTx are passed as an argument, see below:
| # # Name: spams.lasso # # Usage: spams.lasso(X,D= NULL,Q = NULL,q = NULL,return_reg_path = FALSE,L= -1,lambda1= NULL, # lambda2= 0.,mode= 'PENALTY',pos= FALSE,ols= FALSE,numThreads= -1, # max_length_path= -1,verbose=FALSE,cholesky= FALSE) # # Description: # spams.lasso is an efficient implementation of the # homotopy-LARS algorithm for solving the Lasso. # # If the function is called this way spams.lasso(X,D = D, Q = NULL,...), # it aims at addressing the following problems # for all columns x of X, it computes one column alpha of A # that solves # 1) when mode=0 # min_{alpha} ||x-Dalpha||_2^2 s.t. ||alpha||_1 <= lambda1 # 2) when mode=1 # min_{alpha} ||alpha||_1 s.t. ||x-Dalpha||_2^2 <= lambda1 # 3) when mode=2 # min_{alpha} 0.5||x-Dalpha||_2^2 + lambda1||alpha||_1 +0.5 lambda2||alpha||_2^2 # # If the function is called this way spams.lasso(X,D = NULL, Q = Q, q = q,...), # it solves the above optimisation problem, when Q=D'D and q=D'x. # # Possibly, when pos=true, it solves the previous problems # with positivity constraints on the vectors alpha # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # Q: p x p double matrix (Q = D'D) # q: p x n double matrix (q = D'X) # verbose: verbose mode # return_reg_path: # if true the function will return 2 matrices in a list. # lambda1: (parameter) # lambda2: (optional parameter for solving the Elastic-Net) # for mode=0 and mode=1, it adds a ridge on the Gram Matrix # L: (optional), maximum number of steps of the homotopy algorithm (can # be used as a stopping criterion) # pos: (optional, adds non-negativity constraints on the # coefficients, false by default) # mode: (see above, by default: 2) # 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). # cholesky: (optional, default false), choose between Cholesky # implementation or one based on the matrix inversion Lemma # ols: (optional, default false), perform an orthogonal projection # before returning the solution. # max_length_path: (optional) maximum length of the path, by default 4*p # # Output: # A: double sparse p x n matrix (output coefficients) # path: optional, returns the regularisation path for the first signal # A <- spams.lasso(X,return_reg_path = FALSE,...) # v <- spams.lasso(X,return_reg_path = TRUE,...) # A <- v[[1]] # path <- v[[2]] # # Authors: # Julien MAIRAL, 2009 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # # Examples: # m <- 5;n <- 10;nD <- 5 # set.seed(0) # X = matrix(rnorm(m * n),nrow = m,ncol = n,byrow = FALSE) # X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=TRUE) # D = matrix(rnorm(m * nD),nrow = m,ncol = nD,byrow = FALSE) # D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=TRUE) # alpha = spams.lasso(X,D = D,return_reg_path = FALSE,lambda1 = 0.15) # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) .printf("test lasso\n") ############################################## # Decomposition of a large number of signals ############################################## # data generation X = matrix(rnorm(100 * 100000),nrow = 100,ncol = 100000,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) D = matrix(rnorm(100 * 200),nrow = 100,ncol = 200,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) tic = proc.time() alpha = spams.lasso(X,D,return_reg_path = FALSE,lambda1 = 0.15,numThreads = -1,mode = 'PENALTY' ) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) ######################################## # Regularization path of a single signal ######################################## X = matrix(rnorm(64),nrow = 64,ncol = 1,byrow = FALSE) D = matrix(rnorm(64 * 10),nrow = 64,ncol = 10,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) res = spams.lasso(X,D,return_reg_path = TRUE,lambda1 = 0.15,numThreads = -1,mode = 'PENALTY' ) alpha = res[[1]] path = res[[2]] |
This is a fast implementation of a weighted version of LARS [9]. Given a matrix of signals X=[x1,…,xn] in ℝm × n, a matrix of weights W=[w1,…,wn] ∈ ℝp × n, and a dictionary D in ℝm × p, depending on the input parameters, the algorithm returns a matrix of coefficients A=[α1,…,αn] 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
| ||x−Dα||22 s.t. ||diag(w)α||1 ≤ λ, (11) |
or
| ||diag(w)α||1 s.t. ||x−Dα||22 ≤ λ, (12) |
or
|
| ||x−Dα||22 + λ ||diag(w)α||1. (13) |
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-ℓ1 schemes [4].
| # # Name: spams.lassoWeighted. # # Usage: spams.lassoWeighted(X,D,W,L= -1,lambda1= NULL,mode= 'PENALTY',pos= FALSE,numThreads= -1, # verbose=FALSE) # # Description: # spams.lassoWeighted is an efficient implementation of the # LARS algorithm for solving the weighted Lasso. It is optimized # for solving a large number of small or medium-sized # decomposition problem (and not for a single large one). # It first computes the Gram matrix D'D and then perform # a Cholesky-based OMP of the input signals in parallel. # For all columns x of X, and w of W, it computes one column alpha of A # which is the solution of # 1) when mode=0 # min_{alpha} ||x-Dalpha||_2^2 s.t. # ||diag(w)alpha||_1 <= lambda1 # 2) when mode=1 # min_{alpha} ||diag(w)alpha||_1 s.t. # ||x-Dalpha||_2^2 <= lambda1 # 3) when mode=2 # min_{alpha} 0.5||x-Dalpha||_2^2 + # lambda1||diag(w)alpha||_1 # Possibly, when pos=true, it solves the previous problems # with positivity constraints on the vectors alpha # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # W: double p x n matrix (weights) # verbose: verbose mode # lambda1: (parameter) # L: (optional, maximum number of elements of each # decomposition) # pos: (optional, adds positivity constraints on the # coefficients, false by default) # mode: (see above, by default: 2) # 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). # # Output: # A: double sparse p x n matrix (output coefficients) # # Authors: # Julien MAIRAL, 2009 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) .printf("test lasso weighted\n") ############################################## # Decomposition of a large number of signals ############################################## # data generation X = matrix(rnorm(64 * 10000),nrow = 64,ncol = 10000,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) D = matrix(rnorm(64 * 256),nrow = 64,ncol = 256,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) W = matrix(runif(ncol(D) * ncol(X),0,1),nrow = ncol(D),ncol = ncol(X),byrow = FALSE) tic = proc.time() alpha = spams.lassoWeighted(X,D,W,lambda1 = 0.15,numThreads = -1,mode = 'PENALTY') tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) |
This is a variant of spams.lasso with the possibility of adding a mask B=[β1,…,βn], as in spams.ompMask. For every column x of X, β of B, it computes a column α of A, which is the solution of
| ||diag(β)(x−Dα)||22 s.t. ||α||1 ≤ λ, (14) |
or
| ||α||1 s.t. ||diag(β)(x−Dα)||22 ≤ λ |
| , (15) |
or
|
| ||diag(β)(x−Dα)||22 + λ |
| ||α||1 + |
| ||α||22. (16) |
| # # Name: spams.lassoMask # # Usage: spams.lassoMask(X,D,B,L= -1,lambda1= NULL,lambda2= 0.,mode= 'PENALTY',pos= FALSE, # numThreads= -1,verbose=FALSE) # # Description: # spams.lasso is a variant of spams.lasso that handles # binary masks. It aims at addressing the following problems # for all columns x of X, and beta of B, it computes one column alpha of A # that solves # 1) when mode=0 # min_{alpha} ||diag(beta)(x-Dalpha)||_2^2 s.t. ||alpha||_1 <= lambda1 # 2) when mode=1 # min_{alpha} ||alpha||_1 s.t. ||diag(beta)(x-Dalpha)||_2^2 # <= lambda1*||beta||_0/m # 3) when mode=2 # min_{alpha} 0.5||diag(beta)(x-Dalpha)||_2^2 + # lambda1*(||beta||_0/m)*||alpha||_1 + # (lambda2/2)||alpha||_2^2 # Possibly, when pos=true, it solves the previous problems # with positivity constraints on the vectors alpha # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # B: boolean m x n matrix (mask) # p is the number of elements in the dictionary # verbose: verbose mode # lambda1: (parameter) # L: (optional, maximum number of elements of each # decomposition) # pos: (optional, adds positivity constraints on the # coefficients, false by default) # mode: (see above, by default: 2) # lambda2: (optional parameter for solving the Elastic-Net) # for mode=0 and mode=1, it adds a ridge on the Gram Matrix # 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). # # Output: # A: double sparse p x n matrix (output coefficients) # # Authors: # Julien MAIRAL, 2010 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) .printf("test lassoMask\n") ############################################## # Decomposition of a large number of signals ############################################## # data generation X = matrix(rnorm(100 * 100),nrow = 100,ncol = 100,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) D = matrix(rnorm(100 * 20),nrow = 100,ncol = 20,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) mask = (X > 0) # generating a binary mask tic = proc.time() alpha = spams.lassoMask(X,D,mask,lambda1 = 0.15,numThreads = -1,mode = 'PENALTY') tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) |
Coordinate-descent approach for solving Eq. (10) and Eq. (9). Note that unlike spams.lasso, 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.
| # # Name: spams.cd # # Usage: spams.cd(X,D,A0,lambda1 = NULL,mode= 'PENALTY',itermax=100,tol = 0.001,numThreads =-1) # # Description: # spams.cd addresses l1-decomposition problem with a # coordinate descent type of approach. # It is optimized for solving a large number of small or medium-sized # decomposition problem (and not for a single large one). # It first computes the Gram matrix D'D. # This method is particularly well adapted when there is low # correlation between the dictionary elements and when one can benefit # from a warm restart. # It aims at addressing the two following problems # for all columns x of X, it computes a column alpha of A such that # 2) when mode=1 # min_{alpha} ||alpha||_1 s.t. ||x-Dalpha||_2^2 <= lambda1 # For this constraint setting, the method solves a sequence of # penalized problems (corresponding to mode=2) and looks # for the corresponding Lagrange multplier with a simple but # efficient heuristic. # 3) when mode=2 # min_{alpha} 0.5||x-Dalpha||_2^2 + lambda1||alpha||_1 # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # All the columns of D should have unit-norm ! # A0: double sparse p x n matrix (initial guess) # lambda1: (parameter) # mode: (optional, see above, by default 2) # itermax: (maximum number of iterations) # tol: (tolerance parameter) # 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). # # Output: # A: double sparse p x n matrix (output coefficients) # # Authors: # Julien MAIRAL, 2009 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha # is double precision) # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) X = matrix(rnorm(64 * 100),nrow = 64,ncol = 100,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) D = matrix(rnorm(64 * 100),nrow = 64,ncol = 100,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) # parameter of the optimization procedure are chosen lambda1 = 0.015 mode = 'PENALTY' tic = proc.time() alpha = spams.lasso(X,D,lambda1 = lambda1,mode = mode,numThreads = 4) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second for LARS\n", (ncol(X) / t)) E = mean(0.5 * colSums((X - D %*% alpha) ^ 2) + lambda1 * colSums(abs(alpha))) .printf("Objective function for LARS: %g\n",E) tol = 0.001 itermax = 1000 A0 = as(matrix(c(0),nrow = nrow(alpha),ncol = ncol(alpha)),'CsparseMatrix') tic = proc.time() alpha2 = spams.cd(X,D,A0,lambda1 = lambda1,mode = mode,tol = tol, itermax = itermax,numThreads = 4) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second for CD\n", (ncol(X) / t)) E = mean(0.5 * colSums((X - D %*% alpha2) ^ 2) + lambda1 * colSums(abs(alpha2))) .printf("Objective function for CD: %g\n",E) .printf("With Random Design, CD can be much faster than LARS\n") |
This is a fast implementation of the Simultaneous Orthogonal Matching Pursuit algorithm. Given a set of matrices X=[X1,…,Xn] in ℝm × N, where the Xi’s are in ℝm × ni, and a dictionary D in ℝm × p, the algorithm returns a matrix of coefficients A=[A1,…,An] in ℝp × N which is an approximate solution of the following NP-hard problem
| ∀ i |
| ||Xi−DAi||F2 s.t. ||Ai||0,∞ ≤ L. (17) |
or
| ∀ i |
| ||Ai||0,∞ s.t. ||Xi−DAi||F2 ≤ ε ni. (18) |
To be efficient, the method first compute the covariance matrix DTD, then for each signal, it computes DTXi and performs the decomposition with a Cholesky-based algorithm.
| # # Name: spams.somp # (this function has not been intensively tested). # # Usage: spams.somp(X,D,list_groups,L = NULL,eps = 0.,numThreads = -1) # # Description: # spams.somp is an efficient implementation of a # Simultaneous Orthogonal Matching Pursuit algorithm. It is optimized # for solving a large number of small or medium-sized # decomposition problem (and not for a single large one). # It first computes the Gram matrix D'D and then perform # a Cholesky-based OMP of the input signals in parallel. # It aims at addressing the following NP-hard problem # # X is a matrix structured in groups of signals, which we denote # by X=[X_1,...,X_n] # # for all matrices X_i of X, # min_{A_i} ||A_i||_{0,infty} s.t ||X_i-D A_i||_2^2 <= eps*n_i # where n_i is the number of columns of X_i # # or # # min_{A_i} ||X_i-D A_i||_2^2 s.t. ||A_i||_{0,infty} <= L # # Inputs: # X: double m x N matrix (input signals) # m is the signal size # N is the total number of signals # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # All the columns of D should have unit-norm ! # list_groups : int32 vector containing the indices (starting at 0) # of the first elements of each groups. # L: (maximum number of elements in each decomposition) # eps: (threshold on the squared l2-norm of the residual # 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). # # Output: # alpha: double sparse p x N matrix (output coefficients) # # Authors: # Julien MAIRAL, 2010 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) X = matrix(rnorm(64 * 10000),nrow = 64,ncol = 10000,byrow = FALSE) D = matrix(rnorm(64 * 200),nrow = 64,ncol = 200,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) ind_groups = as.vector(seq(from = 0,to = 9999,by = 10),mode= 'integer') tic = proc.time() alpha = spams.somp(X,D,ind_groups,L = 10,eps = 0.1,numThreads=-1) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) |
This is a fast implementation of a simultaneous signal decomposition formulation. Given a set of matrices X=[X1,…,Xn] in ℝm × N, where the Xi’s are in ℝm × ni, and a dictionary D in ℝm × p, the algorithm returns a matrix of coefficients A=[A1,…,An] in ℝp × N which is an approximate solution of the following NP-hard problem
| ∀ i |
| ||Xi−DAi||F2 s.t. ||Ai||1,2 ≤ |
| . (19) |
or
| ∀ i |
| ||Ai||1,2 s.t. ||Xi−DAi||F2 ≤ λ ni. (20) |
To be efficient, the method first compute the covariance matrix DTD, then for each signal, it computes DTXi and performs the decomposition with a Cholesky-based algorithm.
| # # Name: spams.l1L2BCD # (this function has not been intensively tested). # # Usage: spams.l1L2BCD(X,D,alpha0,list_groups,lambda1 = NULL,mode= 'PENALTY',itermax = 100, # tol = 1e-3,numThreads = -1) # # Description: # spams.l1L2BCD is a solver for a # Simultaneous signal decomposition formulation based on block # coordinate descent. # # X is a matrix structured in groups of signals, which we denote # by X=[X_1,...,X_n] # # if mode=2, it solves # for all matrices X_i of X, # min_{A_i} 0.5||X_i-D A_i||_2^2 + lambda1/sqrt(n_i)||A_i||_{1,2} # where n_i is the number of columns of X_i # if mode=1, it solves # min_{A_i} ||A_i||_{1,2} s.t. ||X_i-D A_i||_2^2 <= n_i lambda1 # # Inputs: # X: double m x N matrix (input signals) # m is the signal size # N is the total number of signals # D: double m x p matrix (dictionary) # p is the number of elements in the dictionary # alpha0: double dense p x N matrix (initial solution) # list_groups : int32 vector containing the indices (starting at 0) # of the first elements of each groups. # lambda1: (regularization parameter) # mode: (see above, by default 2) # itermax: (maximum number of iterations, by default 100) # tol: (tolerance parameter, by default 0.001) # 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). # # Output: # alpha: double sparse p x N matrix (output coefficients) # # Authors: # Julien MAIRAL, 2010 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting (even though the output alpha is double # precision) # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) X = matrix(rnorm(64 * 100),nrow = 64,ncol = 100,byrow = FALSE) D = matrix(rnorm(64 * 200),nrow = 64,ncol = 200,byrow = FALSE) D = D / matrix(rep(sqrt(colSums(D*D)),nrow(D)),nrow(D),ncol(D),byrow=T) ind_groups = as.vector(seq(from = 0,to = ncol(X) - 1,by = 10),mode= 'integer')#indices of the first signals in each group itermax = 100 tol = 1e-3 mode = 'PENALTY' lambda1 = 0.15 # squared norm of the residual should be less than 0.1 numThreads = -1 # number of processors/cores to use the default choice is -1 # and uses all the cores of the machine alpha0 = matrix(c(0),nrow = ncol(D), ncol = ncol(X),byrow = FALSE) tic = proc.time() alpha = spams.l1L2BCD(X,D,alpha0,ind_groups,lambda1 = lambda1,mode = mode,itermax = itermax,tol = tol,numThreads = numThreads) tac = proc.time() t = (tac - tic)[['elapsed']] .printf("%f signals processed per second\n",as.double(ncol(X)) / t) |
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=[u1,…,un] in ℝm × n, it finds a matrix V=[v1,…,vn] in ℝm × n so that for all column u of U, it computes a column v of V solving
| ||u−v||22 s.t. ||v||1 ≤ τ, (21) |
or
| ||u−v||22 s.t. λ1||v||1 +λ2||v||22≤ τ, (22) |
or
| ||u−v||22 s.t. λ1||v||1 +λ2||v||22+ λ3 FL(v) ≤ τ, (23) |
or
|
| ||u−v||22 + λ1||v||1 +λ2||v||22+ λ3 FL(v). (24) |
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.
| # # Name: spams.sparseProject # # Usage: spams.sparseProject(U,thrs = 1.0,mode = 1,lambda1 = 0.0,lambda2 = 0.0,lambda3 = 0.0, # pos = FALSE,numThreads = -1) # # Description: # spams.sparseProject solves various optimization # problems, including projections on a few convex sets. # It aims at addressing the following problems # for all columns u of U in parallel # 1) when mode=1 (projection on the l1-ball) # min_v ||u-v||_2^2 s.t. ||v||_1 <= thrs # 2) when mode=2 # min_v ||u-v||_2^2 s.t. ||v||_2^2 + lamuda1||v||_1 <= thrs # 3) when mode=3 # min_v ||u-v||_2^2 s.t ||v||_1 + 0.5lamuda1||v||_2^2 <= thrs # 4) when mode=4 # min_v 0.5||u-v||_2^2 + lamuda1||v||_1 s.t ||v||_2^2 <= thrs # 5) when mode=5 # min_v 0.5||u-v||_2^2 + lamuda1||v||_1 +lamuda2 FL(v) + ... # 0.5lamuda_3 ||v||_2^2 # where FL denotes a "fused lasso" regularization term. # 6) when mode=6 # min_v ||u-v||_2^2 s.t lamuda1||v||_1 +lamuda2 FL(v) + ... # 0.5lamuda3||v||_2^2 <= thrs # # When pos=true and mode <= 4, # it solves the previous problems with positivity constraints # # Inputs: # U: double m x n matrix (input signals) # m is the signal size # n is the number of signals to project # thrs: (parameter) # lambda1: (parameter) # lambda2: (parameter) # lambda3: (parameter) # mode: (see above) # pos: (optional, 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). # # Output: # V: double m x n matrix (output matrix) # # Authors: # Julien MAIRAL, 2009 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # # Note: # this function admits a few experimental usages, which have not # been extensively tested: # - single precision setting # |
The following piece of code illustrates how to use this function. The following piece of code contains usage examples:
| library(spams) set.seed(0) X = matrix(rnorm(20000 * 100),nrow = 20000,ncol = 100,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) .printf( "\n Projection on the l1 ball\n") tic = proc.time() X1 = spams.sparseProject(X,numThreads = -1, pos = FALSE,mode= 1,thrs = 2) tac = proc.time() t = (tac - tic)[['elapsed']] .printf(" Time : %f\n", t) if (t != 0) .printf("%f signals of size %d projected per second\n",(ncol(X) / t),nrow(X)) s = colSums(abs(X1)) .printf("Checking constraint: %f, %f\n",min(s),max(s)) .printf("\n Projection on the Elastic-Net\n") lambda1 = 0.15 tic = proc.time() X1 = spams.sparseProject(X,numThreads = -1, pos = FALSE,mode= 2,thrs = 2,lambda1= lambda1) tac = proc.time() t = (tac - tic)[['elapsed']] .printf(" Time : %f\n", t) if (t != 0) .printf("%f signals of size %d projected per second\n",(ncol(X) / t),nrow(X)) constraints = colSums(X1*X1) + lambda1 * colSums(abs(X1)) .printf("Checking constraint: %f, %f (Projection is approximate : stops at a kink)\n",min(constraints),max(constraints)) .printf("\n Projection on the FLSA\n") lambda1 = 0.7 lambda2 = 0.7 lambda3 = 0.7 X = matrix(runif(2000 * 100,0,1),nrow = 2000,ncol = 100,byrow = FALSE) X = X / matrix(rep(sqrt(colSums(X*X)),nrow(X)),nrow(X),ncol(X),byrow=T) tic = proc.time() X1 = spams.sparseProject(X,numThreads = -1, pos = FALSE,mode= 6,thrs = 2,lambda1= lambda1,lambda2= lambda2,lambda3= lambda3) tac = proc.time() t = (tac - tic)[['elapsed']] .printf(" Time : %f\n", t) if (t != 0) .printf("%f signals of size %d projected per second\n",ncol(X) / t,nrow(X)) # constraints = 0.5 * param[['lambda3']] * colSums(X1*X1) + param[['lambda1']] * colSums(abs(X1)) + param[['lambda2']] * abs(X1[3:nrow(X1),]) - colSums(X1[2,]) constraints = 0.5 * lambda3 * colSums(X1*X1) + lambda1 * colSums(abs(X1)) + lambda2 * abs(X1[3:nrow(X1),]) - colSums(matrix(X1[2,],nrow=1)) .printf("Checking constraint: %f, %f (Projection is approximate : stops at a kink)\n",min(constraints),max(constraints)) |
The function implements an active-set algorithm [37] for solving
|
| ||x−Dα||22 s.t. α ≥ 0 and |
| α[j] = 1. |
| # # Name: spams.decompSimplex # # Usage: spams.decompSimplex(X,Z,computeXtX = FALSE,numThreads = -1) # # Description: # documentation to appear soon # # Inputs: # X: double m x n matrix (input signals) # m is the signal size # n is the number of signals to decompose # # Output: # Z: double % # # Authors: # Yuansi Chen and Julien MAIRAL, 2014 (spams, matlab interface and documentation) # Jean-Paul CHIEZE 2011-2012 (R interface) # |