Previous Up Next

3  Dictionary Learning and Matrix Factorization Toolbox

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

3.1  Function spams.trainDL

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

 
min
D ∈ C
 
 
lim
n → +∞
 
1
n
 
n
∑
i=1
 
 
min
αi
 ⎛
⎜
⎝
1
2
 ||xi−Dαi||22 + ψ(αi)⎞
⎟
⎠		.     (1)

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

#
# Name: trainDL
#
# Usage: spams.trainDL(X,return_model= False,model= None,D = None,numThreads = -1,batchsize = -1,
#               K= -1,lambda1= None,lambda2= 10e-10,iter=-1,t0=1e-5,mode=spams_wrap.PENALTY,
#               posAlpha=False,posD=False,expand=False,modeD=spams_wrap.L2,whiten=False,
#               clean=True,verbose=True,gamma1=0.,gamma2=0.,rho=1.0,iter_updateD=None,
#               stochastic_deprecated=False,modeParam=0,batch=False,log_deprecated=False,
#               logName='')
#
# Description:
#     trainDL is an efficient implementation of the
#     dictionary learning technique presented in
#     
#     "Online Learning for Matrix Factorization and Sparse Coding"
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     arXiv:0908.0050
#     
#     "Online Dictionary Learning for Sparse Coding"      
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     ICML 2009.
#     
#     Note that if you use mode=1 or 2, if the training set has a
#     reasonable size and you have enough memory on your computer, you 
#     should use trainDL_Memory instead.
#     
#     
#     It addresses the dictionary learning problems
#        1) if mode=0
#     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2  s.t. ...
#                                                  ||alpha_i||_1 <= lambda1
#        2) if mode=1
#     min_{D in C} (1/n) sum_{i=1}^n  ||alpha_i||_1  s.t.  ...
#                                           ||x_i-Dalpha_i||_2^2 <= lambda1
#        3) if mode=2
#     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2 + ... 
#                                  lambda1||alpha_i||_1 + lambda1_2||alpha_i||_2^2
#        4) if mode=3, the sparse coding is done with OMP
#     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2  s.t. ... 
#                                                  ||alpha_i||_0 <= lambda1
#        5) if mode=4, the sparse coding is done with OMP
#     min_{D in C} (1/n) sum_{i=1}^n  ||alpha_i||_0  s.t.  ...
#                                           ||x_i-Dalpha_i||_2^2 <= lambda1
#        6) if mode=5, the sparse coding is done with OMP
#     min_{D in C} (1/n) sum_{i=1}^n 0.5||x_i-Dalpha_i||_2^2 +lambda1||alpha_i||_0  
#     
#     
#     C is a convex set verifying
#        1) if modeD=0
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 <= 1 }
#        2) if modeD=1
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... 
#                                                  gamma1||d_j||_1 <= 1 }
#        3) if modeD=2
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... 
#                                  gamma1||d_j||_1 + gamma2 FL(d_j) <= 1 }
#        4) if modeD=3
#           C={  D in Real^{m x p}  s.t.  forall j,  (1-gamma1)||d_j||_2^2 + ... 
#                                  gamma1||d_j||_1 <= 1 }
#                                  
#     Potentially, n can be very large with this algorithm.
#
# Inputs:
#       X:  double m x n matrix   (input signals)
#               m is the signal size
#               n is the number of signals to decompose
#       return_model:     
#               if true the function will return the model
#               as a named list ('A' = A, 'B' = B, 'iter' = n)
#       model:        None or model (as A,B,iter) to use as initialisation
#       D: (optional) double m x p matrix   (dictionary)
#         p is the number of elements in the dictionary
#         When D is not provided, the dictionary is initialized 
#         with random elements from the training set.
#       K: (size of the dictionary, optional is D is provided)
#       lambda1:  (parameter)
#       lambda2:  (optional, by default 0)
#       iter: (number of iterations).  If a negative number is 
#          provided it will perform the computation during the
#          corresponding number of seconds. For instance iter=-5
#          learns the dictionary during 5 seconds.
#       mode: (optional, see above, by default 2) 
#       posAlpha: (optional, adds positivity constraints on the
#         coefficients, false by default, not compatible with 
#         mode =3,4)
#       modeD: (optional, see above, by default 0)
#       posD: (optional, adds positivity constraints on the 
#         dictionary, false by default, not compatible with 
#         modeD=2)
#       gamma1: (optional parameter for modeD >= 1)
#       gamma2: (optional parameter for modeD = 2)
#       batchsize: (optional, size of the minibatch, by default 
#          512)
#       iter_updateD: (optional, number of BCD iterations for the dictionary
#          update step, by default 1)
#       modeParam: (optimization mode).
#          1) if modeParam=0, the optimization uses the 
#             parameter free strategy of the ICML paper
#          2) if modeParam=1, the optimization uses the 
#             parameters rho as in arXiv:0908.0050
#          3) if modeParam=2, the optimization uses exponential 
#             decay weights with updates of the form 
#             A_{t} <- rho A_{t-1} + alpha_t alpha_t^T
#       rho: (optional) tuning parameter (see paper arXiv:0908.0050)
#       t0: (optional) tuning parameter (see paper arXiv:0908.0050)
#       clean: (optional, true by default. prunes 
#          automatically the dictionary from unused elements).
#       verbose: (optional, true by default, increase verbosity)
#       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).
#       expand:    undocumented; modify at your own risks!
#       whiten:    undocumented; modify at your own risks!
#       stochastic_deprecated:    undocumented; modify at your own risks!
#       batch:    undocumented; modify at your own risks!
#       log_deprecated:    undocumented; modify at your own risks!
#       logName:    undocumented; modify at your own risks!
#
# Output:
#       D:     double m x p matrix   (dictionary)
#       model:  the model as A B iter
#        D = spams.trainDL(X,return_model = False,...)
#        (D,model) = spams.trainDL(X,return_model = True,...)
#
# Authors:
# Julien MAIRAL, 2009 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (python 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:

import spams
import numpy as np
img_file = 'boat.png'
try:
    img = Image.open(img_file)
except:
    print("Cannot load image %s : skipping test" %img_file)
I = np.array(img) / 255.
if I.ndim == 3:
    A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])))
    rgb = True
else:
    A = np.asfortranarray(I)
    rgb = False

m = 8
n = 8
X = spams.im2col_sliding(A,m,n,rgb)

X = X - np.tile(np.mean(X,0),(X.shape[0],1))
X = np.asfortranarray(X / np.tile(np.sqrt((X * X).sum(axis=0)),(X.shape[0],1)),dtype = myfloat)
param = { 'K' : 100, # learns a dictionary with 100 elements
          'lambda1' : 0.15, 'numThreads' : 4, 'batchsize' : 400,
          'iter' : 1000}

########## FIRST EXPERIMENT ###########
tic = time.time()
D = spams.trainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

##param['approx'] = 0
# save dictionnary as dict.png
_objective(X,D,param,'dict')

#### SECOND EXPERIMENT ####
print("*********** SECOND EXPERIMENT ***********")

X1 = X[:,0:X.shape[1]//2]
X2 = X[:,X.shape[1]//2 -1:]
param['iter'] = 500
tic = time.time()
(D,model) = spams.trainDL(X1,return_model = True,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f\n' %t)

_objective(X,D,param,'dict1')

# Then reuse the learned model to retrain a few iterations more.
param2 = param.copy()
param2['D'] = D
tic = time.time()
(D,model) = spams.trainDL(X2,return_model = True,model = model,**param2)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)
_objective(X,D,param,'dict2')

#################### THIRD & FOURTH EXPERIMENT ######################
# let us add sparsity to the dictionary itself

print('*********** THIRD EXPERIMENT ***********')
param['modeParam'] = 0
param['iter'] = 1000
param['gamma1'] = 0.3
param['modeD'] = 1

tic = time.time()
D = spams.trainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)
_objective(X,D,param)

print('*********** FOURTH EXPERIMENT ***********')
param['modeParam'] = 0
param['iter'] = 1000
param['gamma1'] = 0.3
param['modeD'] = 3

tic = time.time()
D = spams.trainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)
_objective(X,D,param)

3.2  Function spams.trainDL_Memory

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

#
# Name: trainDL_Memory
#
# Usage: spams.trainDL_Memory(X,D = None,numThreads = -1,batchsize = -1,K= -1,lambda1= None,iter=-1,
#                      t0=1e-5,mode=spams_wrap.PENALTY,posD=False,expand=False,
#                      modeD=spams_wrap.L2,whiten=False,clean=True,gamma1=0.,gamma2=0.,
#                      rho=1.0,iter_updateD=1,stochastic_deprecated=False,modeParam=0,
#                      batch=False,log_deprecated=False,logName='')
#
# Description:
#     trainDL_Memory is an efficient but memory consuming 
#     variant of the dictionary learning technique presented in
#     
#     "Online Learning for Matrix Factorization and Sparse Coding"
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     arXiv:0908.0050
#     
#     "Online Dictionary Learning for Sparse Coding"      
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     ICML 2009.
#     
#     Contrary to the approaches above, the algorithm here 
#        does require to store all the coefficients from all the training
#        signals. For this reason this variant can not be used with large
#        training sets, but is more efficient than the regular online
#        approach for training sets of reasonable size.
#        
#     It addresses the dictionary learning problems
#        1) if mode=1
#     min_{D in C} (1/n) sum_{i=1}^n  ||alpha_i||_1  s.t.  ...
#                                         ||x_i-Dalpha_i||_2^2 <= lambda1
#        2) if mode=2
#     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2 + ... 
#                                                      lambda1||alpha_i||_1  
#                                                      
#     C is a convex set verifying
#        1) if modeD=0
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 <= 1 }
#        1) if modeD=1
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... 
#                                                  gamma1||d_j||_1 <= 1 }
#        1) if modeD=2
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... 
#                                  gamma1||d_j||_1 + gamma2 FL(d_j) <= 1 }
#                                  
#     Potentially, n can be very large with this algorithm.
#
# Inputs:
#       X:  double m x n matrix   (input signals)
#               m is the signal size
#               n is the number of signals to decompose
#       D: (optional) double m x p matrix   (dictionary)
#         p is the number of elements in the dictionary
#         When D is not provided, the dictionary is initialized 
#         with random elements from the training set.
#       K: (size of the dictionary, optional is D is provided)
#       lambda1:  (parameter)
#       iter: (number of iterations).  If a negative number is 
#          provided it will perform the computation during the
#          corresponding number of seconds. For instance iter=-5
#          learns the dictionary during 5 seconds.
#       mode: (optional, see above, by default 2) 
#       modeD: (optional, see above, by default 0)
#       posD: (optional, adds positivity constraints on the 
#         dictionary, false by default, not compatible with 
#         modeD=2)
#       gamma1: (optional parameter for modeD >= 1)
#       gamma2: (optional parameter for modeD = 2)
#       batchsize: (optional, size of the minibatch, by default 
#         512)
#       iter_updateD: (optional, number of BCD iterations for the dictionary 
#           update step, by default 1)
#       modeParam: (optimization mode).
#         1) if modeParam=0, the optimization uses the 
#            parameter free strategy of the ICML paper
#         2) if modeParam=1, the optimization uses the 
#            parameters rho as in arXiv:0908.0050
#         3) if modeParam=2, the optimization uses exponential 
#            decay weights with updates of the form 
#            A_{t} <- rho A_{t-1} + alpha_t alpha_t^T
#       rho: (optional) tuning parameter (see paper arXiv:0908.0050)
#       t0: (optional) tuning parameter (see paper arXiv:0908.0050)
#       clean: (optional, true by default. prunes 
#         automatically the dictionary from unused elements).
#       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).
#       expand:    undocumented; modify at your own risks!
#       whiten:    undocumented; modify at your own risks!
#       stochastic_deprecated:    undocumented; modify at your own risks!
#       batch:    undocumented; modify at your own risks!
#       log_deprecated:    undocumented; modify at your own risks!
#       logName:    undocumented; modify at your own risks!
#
# Output:
#       D:     double m x p matrix   (dictionary)
#       model:  the model as A B iter
#        D = spams.trainDL_Memory(X,...)
#
# Authors:
# Julien MAIRAL, 2009 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (python 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:

import spams
import numpy as np
img_file = 'lena.png'
try:
    img = Image.open(img_file)
except:
    print("Cannot load image %s : skipping test" %img_file)
I = np.array(img) / 255.
if I.ndim == 3:
    A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])))
    rgb = True
else:
    A = np.asfortranarray(I)
    rgb = False

m = 8
n = 8
X = spams.im2col_sliding(A,m,n,rgb)

X = X - np.tile(np.mean(X,0),(X.shape[0],1))
X = np.asfortranarray(X / np.tile(np.sqrt((X * X).sum(axis=0)),(X.shape[0],1)))
X = np.asfortranarray(X[:,np.arange(0,X.shape[1],10)],dtype = myfloat)

param = { 'K' : 200, # learns a dictionary with 100 elements
      'lambda1' : 0.15, 'numThreads' : 4,
      'iter' : 100}

############# FIRST EXPERIMENT  ##################
tic = time.time()
D = spams.trainDL_Memory(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

print('Evaluating cost function...')
lparam = _extract_lasso_param(param)
alpha = spams.lasso(X,D = D,**lparam)
xd = X - D * alpha
R = np.mean(0.5 * (xd * xd).sum(axis=0) + param['lambda1'] * np.abs(alpha).sum(axis=0))
print("objective function: %f" %R)

############# SECOND EXPERIMENT  ##################
tic = time.time()
D = spams.trainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)
print('Evaluating cost function...')
alpha = spams.lasso(X,D = D,**lparam)
xd = X - D * alpha
R = np.mean(0.5 * (xd * xd).sum(axis=0) + param['lambda1'] * np.abs(alpha).sum(axis=0))
print("objective function: %f" %R)

3.3  Function spams.structTrainDL

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

#
# Name: structTrainDL
#
# Usage: spams.structTrainDL(X,return_model= False,model= None,D = None,graph = None,tree = None,
#                     numThreads = -1,tol = 0.000001,fixed_step = True,ista = False,
#                     batchsize = -1,K= -1,lambda1= None,lambda2= 10e-10,lambda3 = 0.,
#                     iter=-1,t0=1e-5,regul = "none",posAlpha=False,posD=False,expand=False,
#                     modeD=spams_wrap.L2,whiten=False,clean=True,verbose=True,gamma1=0.,
#                     gamma2=0.,rho=1.0,iter_updateD=None,stochastic_deprecated=False,
#                     modeParam=0,batch=False,log_deprecated=False,logName='')
#
# Description:
#     structTrainDL is an efficient implementation of the
#     dictionary learning technique presented in
#     
#     "Online Learning for Matrix Factorization and Sparse Coding"
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     arXiv:0908.0050
#     
#     "Online Dictionary Learning for Sparse Coding"      
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     ICML 2009.
#     
#     
#     It addresses the dictionary learning problems
#        min_{D in C} (1/n) sum_{i=1}^n 0.5||x_i-Dalpha_i||_2^2 + lambda1 psi(alpha)
#        where the regularization function psi depends on regul
#        (see proximalFlat for the description of psi,
#         and regul below for allowed values of regul)
#         
#     C is a convex set verifying
#        1) if modeD=0
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 <= 1 }
#        2) if modeD=1
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... 
#                                                  gamma1||d_j||_1 <= 1 }
#        3) if modeD=2
#           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... 
#                                  gamma1||d_j||_1 + gamma2 FL(d_j) <= 1 }
#        4) if modeD=3
#           C={  D in Real^{m x p}  s.t.  forall j,  (1-gamma1)||d_j||_2^2 + ... 
#                                  gamma1||d_j||_1 <= 1 }
#                                  
#     Potentially, n can be very large with this algorithm.
#
# Inputs:
#       X:  double m x n matrix   (input signals)
#               m is the signal size
#               n is the number of signals to decompose
#       return_model:     
#               if true the function will return the model
#               as a named list ('A' = A, 'B' = B, 'iter' = n)
#       model:        None or model (as A,B,iter) to use as initialisation
#       D: (optional) double m x p matrix   (dictionary)
#         p is the number of elements in the dictionary
#         When D is not provided, the dictionary is initialized 
#         with random elements from the training set.
#       K: (size of the dictionary, optional is D is provided)
#       lambda1:  (parameter)
#       lambda2:  (optional, by default 0)
#       lambda3: (optional, regularization parameter, 0 by default)
#       iter: (number of iterations).  If a negative number is 
#          provided it will perform the computation during the
#          corresponding number of seconds. For instance iter=-5
#          learns the dictionary during 5 seconds.
#       regul: choice of regularization : one of
#           'l0' 'l1' 'l2' 'linf' 'none' 'elastic-net' 'fused-lasso'
#           'graph' 'graph-ridge' 'graph-l2' 'tree-l0' 'tree-l2' 'tree-linf' 
#       tree: struct (see documentation of proximalTree);
#           needed for regul of graph kind.
#       graph: struct (see documentation of proximalGraph);
#           needed for regul of tree kind.
#       posAlpha: (optional, adds positivity constraints on the
#           coefficients, false by default.
#       modeD: (optional, see above, by default 0)
#       posD: (optional, adds positivity constraints on the 
#         dictionary, false by default, not compatible with 
#         modeD=2)
#       gamma1: (optional parameter for modeD >= 1)
#       gamma2: (optional parameter for modeD = 2)
#       batchsize: (optional, size of the minibatch, by default 
#          512)
#       iter_updateD: (optional, number of BCD iterations for the dictionary
#          update step, by default 1)
#       modeParam: (optimization mode).
#          1) if modeParam=0, the optimization uses the 
#             parameter free strategy of the ICML paper
#          2) if modeParam=1, the optimization uses the 
#             parameters rho as in arXiv:0908.0050
#          3) if modeParam=2, the optimization uses exponential 
#             decay weights with updates of the form 
#             A_{t} <- rho A_{t-1} + alpha_t alpha_t^T
#       ista: (optional, use ista instead of fista, false by default).
#       tol: (optional, tolerance for stopping criteration, which is a relative duality gap
#       fixed_step: (deactive the line search for L in fista and use K instead)
#       rho: (optional) tuning parameter (see paper arXiv:0908.0050)
#       t0: (optional) tuning parameter (see paper arXiv:0908.0050)
#       clean: (optional, true by default. prunes 
#          automatically the dictionary from unused elements).
#       verbose: (optional, true by default, increase verbosity)
#       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).
#       expand:    undocumented; modify at your own risks!
#       whiten:    undocumented; modify at your own risks!
#       stochastic_deprecated:    undocumented; modify at your own risks!
#       batch:    undocumented; modify at your own risks!
#       log_deprecated:    undocumented; modify at your own risks!
#       logName:    undocumented; modify at your own risks!
#
# Output:
#       D:     double m x p matrix   (dictionary)
#       model:  the model as A B iter
#        D = spams.structTrainDL(X,return_model = False,...)
#        (D,model) = spams.structTrainDL(X,return_model = True,...)
#
# Authors:
# Julien MAIRAL, 2009 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (python 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:

import spams
import numpy as np
img_file = 'lena.png'
try:
    img = Image.open(img_file)
except Exception as e:
    print("Cannot load image %s (%s) : skipping test" %(img_file,e))
I = np.array(img) / 255.
if I.ndim == 3:
    A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])),dtype = myfloat)
    rgb = True
else:
    A = np.asfortranarray(I,dtype = myfloat)
    rgb = False

m = 8
n = 8
X = spams.im2col_sliding(A,m,n,rgb)

X = X - np.tile(np.mean(X,0),(X.shape[0],1))
X = np.asfortranarray(X / np.tile(np.sqrt((X * X).sum(axis=0)),(X.shape[0],1)),dtype = myfloat)
param = { 'K' : 64, # learns a dictionary with 100 elements
          'lambda1' : 0.05, 'tol' : 1e-3,
          'numThreads' : 4, 'batchsize' : 400,
          'iter' : 20}
paramL = {'lambda1' : 0.05, 'numThreads' : 4}

param['regul'] = 'l1'
print("with Fista Regression %s" %param['regul'])
tic = time.time()
D = spams.structTrainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)
_objective(X,D,param)

#
param['regul'] = 'l2'
print("with Fista Regression %s" %param['regul'])
tic = time.time()
D = spams.structTrainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

_objective(X,D,param)

#
param['regul'] = 'elastic-net'
print("with Fista %s" %param['regul'])
param['lambda2'] = 0.1
tic = time.time()
D = spams.structTrainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

_objective(X,D,param)
## if we want a pause :
##    s = raw_input("graph> ")
########### GRAPH
param['lambda1'] = 0.1
param['tol'] = 1e-5
param['K'] = 10
eta_g = np.array([1, 1, 1, 1, 1],dtype=myfloat)

groups = ssp.csc_matrix(np.array([[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]],dtype=np.bool),dtype=np.bool)

groups_var = ssp.csc_matrix(np.array([[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]],dtype=np.bool),dtype=np.bool)

graph = {'eta_g': eta_g,'groups' : groups,'groups_var' : groups_var}
param['graph'] = graph
param['tree'] = None

param['regul'] = 'graph'
print("with Fista %s" %param['regul'])

tic = time.time()
D = spams.structTrainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

_objective(X,D,param)

param['regul'] = 'graph-ridge'
print("with Fista %s" %param['regul'])
tic = time.time()
D = spams.structTrainDL(X,**param)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

_objective(X,D,param)
## if we want a pause :
##    s = raw_input("tree> ")
##### TREE
tree_data = """0 1. [] -> 1 4

3.4  Function nmf

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

#
# Name: nmf
#
# Usage: spams.nmf(X,return_lasso= False,model= None,numThreads = -1,batchsize = -1,K= -1,iter=-1,
#           t0=1e-5,clean=True,rho=1.0,modeParam=0,batch=False)
#
# Description:
#     trainDL is an efficient implementation of the
#     non-negative matrix factorization technique presented in 
#     
#     "Online Learning for Matrix Factorization and Sparse Coding"
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     arXiv:0908.0050
#     
#     "Online Dictionary Learning for Sparse Coding"      
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     ICML 2009.
#     
#     Potentially, n can be very large with this algorithm.
#
# Inputs:
#       X:  double m x n matrix   (input signals)
#               m is the signal size
#               n is the number of signals to decompose
#       return_lasso:     
#                          if true the function will return a tuple of matrices.
#       K: (number of required factors)
#       iter: (number of iterations).  If a negative number 
#         is provided it will perform the computation during the
#         corresponding number of seconds. For instance iter=-5
#         learns the dictionary during 5 seconds.
#       batchsize: (optional, size of the minibatch, by default 
#          512)
#       modeParam: (optimization mode).
#          1) if modeParam=0, the optimization uses the 
#             parameter free strategy of the ICML paper
#          2) if modeParam=1, the optimization uses the 
#             parameters rho as in arXiv:0908.0050
#          3) if modeParam=2, the optimization uses exponential 
#             decay weights with updates of the form  
#             A_{t} <- rho A_{t-1} + alpha_t alpha_t^T
#       rho: (optional) tuning parameter (see paper 
#       arXiv:0908.0050)
#       t0: (optional) tuning parameter (see paper 
#       arXiv:0908.0050)
#       clean: (optional, true by default. prunes automatically 
#         the dictionary from unused elements).
#       batch: (optional, false by default, use batch learning 
#         instead of online learning)
#       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).
#       model: struct (optional) learned model for "retraining" the data.
#
# Output:
#       U: double m x p matrix   
#       V: double p x n matrix   (optional)
#       model: struct (optional) learned model to be used for 
#         "retraining" the data.
#         U = spams.nmf(X,return_lasso = False,...)
#     (U,V) = spams.nmf(X,return_lasso = True,...)
#
# Authors:
# Julien MAIRAL, 2009 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (python interface)
#

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

import spams
import numpy as np
img_file = 'boat.png'
try:
    img = Image.open(img_file)
except:
    print("Cannot load image %s : skipping test" %img_file)
I = np.array(img) / 255.
if I.ndim == 3:
    A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])),dtype = myfloat)
    rgb = True
else:
    A = np.asfortranarray(I,dtype = myfloat)
    rgb = False

m = 16
n = 16
X = spams.im2col_sliding(A,m,n,rgb)
X = X[:,::10]
X = np.asfortranarray(X / np.tile(np.sqrt((X * X).sum(axis=0)),(X.shape[0],1)),dtype = myfloat)
########## FIRST EXPERIMENT ###########
tic = time.time()
(U,V) = spams.nmf(X,return_lasso= True,K = 49,numThreads=4,iter = -5)
tac = time.time()
t = tac - tic
print('time of computation for Dictionary Learning: %f' %t)

print('Evaluating cost function...')
Y = X - U * V
R = np.mean(0.5 * (Y * Y).sum(axis=0))
print('objective function: %f' %R)


# Archetypal Analysis, run first steps with FISTA and run last steps with activeSet,

3.5  Function nnsc

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

#
# Name: nmf
#
# Usage: spams.nnsc(X,return_lasso= False,model= None,lambda1= None,numThreads = -1,batchsize = -1,
#            K= -1,iter=-1,t0=1e-5,clean=True,rho=1.0,modeParam=0,batch=False)
#
# Description:
#     trainDL is an efficient implementation of the
#     non-negative sparse coding technique presented in 
#     
#     "Online Learning for Matrix Factorization and Sparse Coding"
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     arXiv:0908.0050
#     
#     "Online Dictionary Learning for Sparse Coding"      
#     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro
#     ICML 2009.
#     
#     Potentially, n can be very large with this algorithm.
#
# Inputs:
#       X:  double m x n matrix   (input signals)
#               m is the signal size
#               n is the number of signals to decompose
#       return_lasso:     
#                          if true the function will return a tuple of matrices.
#       K: (number of required factors)
#       lambda1: (parameter)
#       iter: (number of iterations).  If a negative number 
#          is provided it will perform the computation during the
#          corresponding number of seconds. For instance iter=-5
#          learns the dictionary during 5 seconds.
#       batchsize: (optional, size of the minibatch, by default 
#          512)
#       modeParam: (optimization mode).
#          1) if modeParam=0, the optimization uses the 
#             parameter free strategy of the ICML paper
#          2) if modeParam=1, the optimization uses the 
#             parameters rho as in arXiv:0908.0050
#          3) if modeParam=2, the optimization uses exponential 
#             decay weights with updates of the form 
#             A_{t} <- rho A_{t-1} + alpha_t alpha_t^T
#       rho: (optional) tuning parameter (see paper
#       arXiv:0908.0050)
#       t0: (optional) tuning parameter (see paper 
#       arXiv:0908.0050)
#       clean: (optional, true by default. prunes automatically 
#          the dictionary from unused elements).
#       batch: (optional, false by default, use batch learning 
#          instead of online learning)
#       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).
#       model: struct (optional) learned model for "retraining" the data.
#
# Output:
#       U: double m x p matrix   
#       V: double p x n matrix   (optional)
#       model: struct (optional) learned model to be used for 
#          "retraining" the data.
#          U = spams.nnsc(X,return_lasso = False,...)
#      (U,V) = spams.nnsc(X,return_lasso = True,...)
#
# Authors:
# Julien MAIRAL, 2009 (spams, matlab interface and documentation)
# Jean-Paul CHIEZE 2011-2012 (python interface)
#

3.6  Function spams.archetypalAnalysis

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

#
# Name: archetypalAnalysis
#
# Usage: spams.archetypalAnalysis(X,p = 10,Z0 = None,returnAB = False,robust=False,epsilon=1e-3,
#                          computeXtX=False,stepsFISTA=3,stepsAS=50,randominit=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 (python interface)
#

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

Previous Up Next