-  3.39.25
pyfaust.fact Namespace Reference

The pyfaust factorization module. More...

Functions

def svdtj (M, nGivens=None, tol=0, err_period=100, relerr=True, nGivens_per_fac=None, enable_large_Faust=False, **kwargs)
 Performs a singular value decomposition and returns the left and right singular vectors as Faust transforms. More...
 
def pinvtj (M, nGivens=None, tol=0, err_period=100, relerr=True, nGivens_per_fac=None, enable_large_Faust=False, **kwargs)
 Computes the pseudoinverse of M using svdtj. More...
 
def eigtj (M, nGivens=None, tol=0, err_period=100, order='ascend', relerr=True, nGivens_per_fac=None, verbosity=0, enable_large_Faust=False)
 Performs an approximate eigendecomposition of M and returns the eigenvalues in W along with the corresponding normalized right eigenvectors (as the columns of the Faust object V). More...
 
def palm4msa (M, p, ret_lambda=False, backend=2016, on_gpu=False)
 Factorizes the matrix M with Palm4MSA algorithm using the parameters set in p. More...
 
def palm4msa_mhtp (M, palm4msa_p, mhtp_p, ret_lambda=False, on_gpu=False)
 Runs the MHTP-PALM4MSA algorithm to factorize the matrix M. More...
 
def hierarchical_mhtp (M, hierar_p, mhtp_p, ret_lambda=False, ret_params=False, on_gpu=False)
 Runs the MHTP-PALM4MSA hierarchical factorization algorithm on the matrix M. More...
 
def hierarchical (M, p, ret_lambda=False, ret_params=False, backend=2016, on_gpu=False)
 Factorizes the matrix M using the hierarchical PALM4MSA algorithm according to the parameters set in p. More...
 
def butterfly (M, type="bbtree", perm=None, diag_opt=False, mul_perm=None)
 Factorizes M according to a butterfly support and optionally a permutation using the algorithms described in [1]. More...
 

Detailed Description

The pyfaust factorization module.

This module gives access to the main factorization algorithms of FAuST. These algorithms can factorize a dense matrix into a sparse product (i.e. a Faust object).

There are several factorization algorithms.

  • The first one is PALM4MSA: which stands for Proximal Alternating Linearized Minimization for Multi-layer Sparse Approximation. Note that Palm4MSA is not intended to be used directly. You should rather rely on the second algorithm.
  • The second one is the Hierarchical Factorization algorithm: this is the central algorithm to factorize a dense matrix into a Faust. It makes iterative use of Palm4MSA to proceed with the factorization of a given dense matrix.
  • The third group of algorithms is for approximate eigenvalue decomposition (eigtj) and singular value decomposition (svdtj).
  • The fourth algorithm is pyfaust.fact.butterfly.

Function Documentation

◆ butterfly()

def pyfaust.fact.butterfly (   M,
  type = "bbtree",
  perm = None,
  diag_opt = False,
  mul_perm = None 
)

Factorizes M according to a butterfly support and optionally a permutation using the algorithms described in [1].

The result is a Faust F of the form BP where B has a butterfly structure and P is a permutation matrix determined by the optional parameter perm.

Parameters
M(numpy ndarray) The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance). M must be square and its dimension must be a power of two.
type(str) the type of factorization 'right'ward, 'left'ward or 'bbtree'. More precisely: if 'left' (resp. 'right') is used then at each stage of the factorization the most left factor (resp. the most right factor) is split in two. If 'bbtree' is used then the matrix is factorized according to a Balanced Binary Tree (which is faster as it allows parallelization).
permfive kinds of values are possible for this argument.
  1. perm is a list of column indices of the permutation matrix P which is such that the returned Faust is F = B@P where B is the Faust butterfly approximation of M @ P.T. If the list of indices is not a valid permutation the behaviour is undefined (however an invalid size or an out of bound index raise an exception).
  2. perm is a list of lists of permutation column indices as defined in 1. In that case, all permutations passed to the function are used as explained in 1, each one producing a Faust, the best one (that is the best approximation of M in the Frobenius norm) is kept and returned by butterfly.
  3. perm is 'default_8', this is a particular case of 2. Eight default permutations are used. For the definition of those permutations please refer to [2].
  4. perm is 'bitrev': in that case the permutation is the bit-reversal permutation (cf. pyfaust.bitrev_perm).
  5. By default this argument is None, no permutation is used (this is equivalent to using the identity permutation matrix in 1).
diag_opt(bool) if True then the returned Faust is optimized using pyfaust.opt_butterfly_faust.
mul_perm(bool) decides if the permutation is multiplied into the rightest butterfly factor (mul_perm=True) or if this permutation is left apart as the rightest factor of the Faust (mul_perm=False). It can't be True if diag_opt is True (an error is raised otherwise). Defaultly, mul_perm=None which implies that mul_perm is True iff diag_opt is False.
Note
Below is an example of how to create a permutation scipy CSR matrix from a permutation list of indices (as defined by the perm argument) and conversely how to convert a permutation matrix to a list of indices of permutation.
>>> from scipy.sparse import random, csr_matrix
>>> from numpy.random import permutation
>>> import numpy as np
>>> np.random.seed(42)
>>> I = permutation(8) # random permutation as a list of indices
>>> I
array([1, 5, 0, 7, 2, 4, 3, 6])
>>> n = len(I)
>>> # convert a permutation as indices to a csr_matrix
>>> P = csr_matrix((np.ones(n), (I, np.arange(n))))
>>> P.toarray()
array([[0., 0., 1., 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., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 1., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1.],
[0., 0., 0., 1., 0., 0., 0., 0.]])
>>> # convert a permutation as a list of indices to a permutation matrix P as a csr_matrix
>>> I_ = P.T.nonzero()[1]
>>> I_
array([1, 5, 0, 7, 2, 4, 3, 6], dtype=int32)
>>> np.allclose(I_, I)
True
Returns
The Faust F which is an approximation of M according to a butterfly support.

Examples

>>> import numpy as np
>>> from random import randint
>>> from pyfaust.fact import butterfly
>>> from pyfaust import Faust, wht, dft
>>> H = wht(8).toarray()
>>> F = butterfly(H, type='bbtree')
>>> # compute the error
>>> (F-H).norm()/Faust(H).norm()
1.2560739454502295e-15
>>> # the same can be done with dft in place of wht
>>> # all you need is to specify the bit-reversal permutation
>>> # since the Discrete Fourier Transform is the product of a butterfly factors with this particular permutation
>>> DFT = dft(8).toarray()
>>> F = butterfly(DFT, type='bbtree', perm='bitrev')
>>> # compute the error
>>> (F-DFT).norm()/Faust(DFT).norm()
1.1427230601405052e-15

Use simple permutations:

>>> import numpy as np
>>> from random import randint
>>> from pyfaust.fact import butterfly
>>> M = np.random.rand(4, 4)
>>> # without any permutation
>>> F1 = butterfly(M, type='bbtree')
>>> # which is equivalent to using the identity permutation
>>> p = np.arange(0, 4)
>>> F2 = butterfly(M, type='bbtree', perm=p)
>>> # compute the relative diff
>>> (F2-F1).norm()/F1.norm()
0.0
>>> # then try another permutation
>>> p2 = [1, 0, 3, 2]
>>> F3 = butterfly(M, type='bbtree', perm=p2)

Use butterfly with a permutation defined by a list of indices J:

>>> import numpy as np
>>> from pyfaust.fact import butterfly
>>> from pyfaust import Faust, wht, dft
>>> H = wht(8).toarray()
>>> J = np.arange(7, -1, -1)
>>> F = butterfly(H, type='bbtree', perm=J)
>>> # this is equivalent to passing a list containing a single permutation :
>>> # F = butterfly(H, type='bbtree', perm=[J])
use butterfly with successive permutations J1 and J2
and keep the best approximation
>>> J1 = J
>>> from numpy.random import permutation
>>> J2 = permutation(len(J)) # another random permutation
>>> F = butterfly(H, type='bbtree', perm=[J1, J2])
>>> # or to to use the 8 default permutations (keeping the best approximation resulting Faust)
>>> F = butterfly(H, type='bbtree', perm='default_8')
References
[1] Quoc-Tung Le, Léon Zheng, Elisa Riccietti, Rémi Gribonval. Fast learning of fast transforms, with guarantees. ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing, May 2022, Singapore, Singapore. (hal-03438881)
[2] T. Dao, A. Gu, M. Eichhorn, A. Rudra, and C. Re, “Learning Fast Algorithms for Linear Transforms Using Butterfly Factorizations,” in Proceedings of the 36th International Conference on Machine Learning. June 2019, pp. 1517–1527, PMLR
See also
pyfaust.wht, pyfaust.dft,
pyfaust.rand_butterfly

◆ eigtj()

def pyfaust.fact.eigtj (   M,
  nGivens = None,
  tol = 0,
  err_period = 100,
  order = 'ascend',
  relerr = True,
  nGivens_per_fac = None,
  verbosity = 0,
  enable_large_Faust = False 
)

Performs an approximate eigendecomposition of M and returns the eigenvalues in W along with the corresponding normalized right eigenvectors (as the columns of the Faust object V).

The output is such that V*numpy.diag(W)*V.H approximates M. V is a product of Givens rotations obtained by truncating the Jacobi algorithm.

The trade-off between accuracy and complexity of V can be set through the parameters nGivens and tol that define the targeted number of Givens rotations and targeted error.

Parameters
M(numpy.ndarray or csr_matrix) the matrix to diagonalize. Must be real and symmetric, or complex hermitian. Can be in dense or sparse format. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
nGivens(int) targeted number of Givens rotations (this argument is optional only if tol is set).
tol(float) the tolerance error at which the algorithm stops. The default value is zero so that stopping is based on reaching the targeted nGivens (this argument is optional only if nGivens is set).
err_period(int) it defines the period, in number of factors of V the error is compared to tol (reducing the period spares some factors but increases slightly the computational cost because the error is computed more often).
order(str) order of eigenvalues, possible choices are ‘ascend, 'descend' or 'undef' (to avoid a sorting operation and save some time).
nGivens_per_fac(int) targeted number of Givens rotations per factor of V. Must be an integer between 1 to floor(M.shape[0]/2) (the default value).
relErr(bool) the type of error used as stopping criterion. True for the relative error norm(V*D*V'-M, 'fro')/norm(M, 'fro'), False for the absolute error norm(V*D*V'-M, 'fro').
verbosity(int) the level of verbosity. The greater the value the more info is displayed. It can be helpful to understand for example why the algorithm stopped before reaching the tol error or the number of Givens (nGivens).
enable_large_Faust(bool) if true, it allows to compute a transform that doesn't worth it regarding its complexity compared to the matrix M. Otherwise by default, an exception is raised before the algorithm starts.
Remarks
  • When ‘nGivens’ and ‘tol’ are used simultaneously, the number of Givens rotations in V may be smaller than specified by ‘nGivens’ if the error criterion is met first, and the achieved error may be larger than specified if ‘nGivens’ is reached first during the iterations of the truncated Jacobi algorithm.
  • When nGivens_per_fac > 1, all factors have exactly nGivens_per_fac except the leftmost one which may have fewer if the total number of Givens rotations is not a multiple of nGivens_per_fact
Returns
  • W: (numpy.ndarray) the vector of the approximate eigenvalues of M (in ascending order by default).
  • V: (Faust) the Faust object representing the approximate eigenvector transform. The column V[:, i] is the eigenvector corresponding to the eigenvalue W[i].
References
[1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast graph Fourier transforms via multi-layer sparse approximations", IEEE Transactions on Signal and Information Processing over Networks 2018, 4(2), pp 407-420 https://hal.inria.fr/hal-01416110

Examples

>>> import numpy as np
>>> from pyfaust.fact import eigtj
>>> from scipy.io import loadmat
>>> from os.path import sep
>>> from pyfaust.demo import get_data_dirpath
>>> from numpy.linalg import norm
>>> # get a graph Laplacian to diagonalize
>>> demo_path = sep.join((get_data_dirpath(),'Laplacian_256_community.mat'))
>>> data_dict = loadmat(demo_path)
>>> Lap = data_dict['Lap'].astype('float64')
>>> Dhat, Uhat = eigtj(Lap, nGivens=Lap.shape[0]*100, enable_large_Faust=True)
>>> # Uhat is the Fourier matrix/eigenvectors approximation as a Faust
>>> # (200 factors)
>>> # Dhat the eigenvalues diagonal matrix approx.
>>> print("err: ", np.round(norm(Lap-Uhat@np.diag(Dhat)@Uhat.H)/norm(Lap), decimals=4)) # about 6.5e-3
err: 0.0065
>>> Dhat2, Uhat2 = eigtj(Lap, tol=0.01)
>>> assert(norm(Lap-Uhat2@np.diag(Dhat2)@Uhat2.H)/norm(Lap) < .011)
>>> # and then asking for an absolute error
>>> Dhat3, Uhat3 = eigtj(Lap, tol=0.1, relerr=False)
>>> assert(norm(Lap-Uhat3@np.diag(Dhat3)@Uhat3.H) < .11)
>>> # now recompute Uhat2, Dhat2 but asking a descending order of eigenvalues
>>> Dhat4, Uhat4 = eigtj(Lap, tol=0.01, order='descend')
>>> assert((Dhat4[::-1] == Dhat2[::]).all())
>>> # and now with no sort
>>> Dhat5, Uhat5 = eigtj(Lap, tol=0.01, order='undef')
>>> assert((np.sort(Dhat5) == Dhat2).all())
See also
svdtj

◆ hierarchical()

def pyfaust.fact.hierarchical (   M,
  p,
  ret_lambda = False,
  ret_params = False,
  backend = 2016,
  on_gpu = False 
)

Factorizes the matrix M using the hierarchical PALM4MSA algorithm according to the parameters set in p.

Note
pyfaust.faust_fact() is an alias of this function.

Basically, the hierarchical factorization works in a sequence of splits in two of the last factor. The first split consists to split the matrix M in two, when p.is_fact_side_left is False (which is the case by default), the factorization works toward right, so M is factorized as follows:

\[M \approx S_1 R_1\]

We call \(S_1\) the main factor and \(R_1\) the residual factor. On step 2, \(R_1\) is split in two such as \(R_1 \approx S_2 R_2\), which gives:

\[M \approx S_1 S_2 R_2 \]

And so on until the algorithm reaches the targeted number of factors:

\[M \approx S_1 S_2 ... S_{N-1} R_N \]

If p.is_fact_side_left is False, the residual is factorized toward left, so it gives rather :

\[M \approx R_1 S_1 \\ \\ M \approx R_2 S_2 S_1 \\ \vdots \\ M \approx R_N S_{N-1} ... S_2 S_1 \]

Parameters
M(numpy array) the array to factorize. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
p(pyfaust.factparams.ParamsHierarchical, list or str) is a set of hierarchical factorization parameters. It might be a fully defined instance of parameters (pyfaust.factparams.ParamsHierarchical) or a simplified expression which designates a pre-defined parameterization:
  • 'hadamard' to use pre-defined parameters typically used to factorize a Hadamard matrix of order a power of two (see pyfaust.demo.hadamard).
  • ['rectmat', j, k, s] to use pre-defined parameters used for instance in factorization of the MEG matrix which is a rectangular matrix of size m*n such that m < n (see pyfaust.demo.bsl); j is the number of factors, k the sparsity of the main factor's columns, and s the sparsity of rows for all other factors except the residual factor (that is the first factor here because the factorization is made toward the left – is_side_fact_left == True, cf. pyfaust.factparams.ParamsHierarchical and pyfaust.factparams.ParamsHierarchicalRectMat).
    The residual factor has a sparsity of P*rho^(num_facts-1).
    By default, rho == .8 and P = 1.4. It's possible to set custom values with for example p == ( ['rectmat', j, k, s], {'rho':.4, 'P':.7 }).
    The sparsity is here the number of nonzero elements.
backend(int) the C++ implementation to use (default to 2016, 2020 backend should be faster for most of the factorizations).
on_gpu(bool) if True the GPU implementation is executed (this option applies only to 2020 backend).
ret_lambda(bool) set to True to ask the function to return the scale factor (False by default).
ret_params(bool) set to True to ask the function to return the ParamsHierarchical instance used (False by default). It is useful for consulting what precisely means the simplified parameterizations used to generate a pyfaust.factparams.ParamsHierarchical instance and possibly adjust its attributes to factorize again.
Note
If backend==2020 and regardless to the StoppingCriterion-s defined in p, it is possible to stop any internal call to PALM4MSA manually at any iteration by the key combination CTRL-C. The last Faust computed in the PALM4MSA instance will be used to continue the hierarchical factorization. A typical use case is when the verbose mode is enabled and you see that the error doesn't change anymore or only slightly, you might stop iterations by typing CTRL-C.
Returns
  • F the Faust object result of the factorization: Faust \(([S_1, S_2, ... ,S_{N-1}, R_N]])\) if p.is_fact_side_left == False, Faust \(([R_N, S_{N-1}, ... , S_2, S_1])\) otherwise.
  • if ret_lambda == True (and ret_params == False), then the function returns the tuple (F,_lambda) (_lambda is the scale factor at the end of factorization).
  • if ret_params == True (and ret_lambda == False), then the function returns the tuple (F, p) (p being the ParamsHierarchical instance really used by the algorithm).
  • if ret_lambda == True and ret_params == True, then the function returns the tuple (F, _lambda, p).


Examples

1. Fully Defined Parameters for a Random Matrix Factorization

>>> from pyfaust.fact import hierarchical
>>> from pyfaust.factparams import ParamsHierarchical, ConstraintList, StoppingCriterion
>>> import numpy as np
>>> M = np.random.rand(500, 32)
>>> fact_cons = ConstraintList('splin', 5, 500, 32, 'sp', 96, 32, 32, 'sp', 96, 32, 32)
>>> res_cons = ConstraintList('normcol', 1, 32, 32, 'sp', 666, 32, 32, 'sp', 333, 32, 32)
>>> # or alternatively using pyfaust.proj
>>> # from pyfaust.proj import *
>>> # res_cons = [normcol((32,32), 1), sp((32,32), 666), sp((32,32), 333)]
>>> stop_crit1 = StoppingCriterion(num_its=200)
>>> stop_crit2 = StoppingCriterion(num_its=200)
>>> param = ParamsHierarchical(fact_cons, res_cons, stop_crit1, stop_crit2)
>>> F = hierarchical(M, param)
>>> F

Faust size 500x32, density 0.189063, nnz_sum 3025, 4 factor(s):

  • FACTOR 0 (double) SPARSE, size 500x32, density 0.15625, nnz 2500
  • FACTOR 1 (double) SPARSE, size 32x32, density 0.09375, nnz 96
  • FACTOR 2 (double) SPARSE, size 32x32, density 0.09375, nnz 96
  • FACTOR 3 (double) SPARSE, size 32x32, density 0.325195, nnz 333

2. Simplified Parameters for Hadamard Factorization

>>> from pyfaust import wht
>>> from pyfaust.fact import hierarchical
>>> from numpy.linalg import norm
>>> # generate a Hadamard Faust of size 32x32
>>> FH = wht(32)
>>> H = FH.toarray() # the full matrix version
>>> # factorize it
>>> FH2 = hierarchical(H, 'hadamard');
>>> # test the relative error
>>> (FH-FH2).norm('fro')/FH.norm('fro') # the result is about 1e-16, the factorization is accurate
1.6883057247219393e-16
>>> FH

Faust size 32x32, density 0.3125, nnz_sum 320, 5 factor(s):

  • FACTOR 0 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 1 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 2 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 3 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 4 (double) SPARSE, size 32x32, density 0.0625, nnz 64
>>> FH2

Faust size 32x32, density 0.3125, nnz_sum 320, 5 factor(s):

  • FACTOR 0 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 1 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 2 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 3 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 4 (double) SPARSE, size 32x32, density 0.0625, nnz 64


See also
pyfaust.factparams.ParamsHierarchicalWHT

3. Simplified Parameters for a Rectangular Matrix Factorization (the BSL demo MEG matrix)

>>> from pyfaust import *
>>> from pyfaust.fact import hierarchical
>>> from scipy.io import loadmat
>>> from pyfaust.demo import get_data_dirpath
>>> import numpy as np
>>> d = loadmat(get_data_dirpath()+'/matrix_MEG.mat')
>>> MEG = d['matrix'].T
>>> num_facts = 9
>>> k = 10
>>> s = 8
>>> MEG16 = hierarchical(MEG, ['rectmat', num_facts, k, s])
>>> MEG16

Faust size 204x8193, density 0.0631655, nnz_sum 105573, 9 factor(s):

  • FACTOR 0 (double) SPARSE, size 204x204, density 0.293613, nnz 12219
  • FACTOR 1 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 2 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 3 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 4 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 5 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 6 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 7 (double) SPARSE, size 204x204, density 0.0392157, nnz 1632
  • FACTOR 8 (double) SPARSE, size 204x8193, density 0.0490196, nnz 81930
>>> # verify the constraint k == 10, on column 4
>>> np.count_nonzero(MEG16.factors(8)[:,4].toarray())
10
>>> # now verify the s constraint is respected on MEG16 factor 1
>>> np.count_nonzero(MEG16.factors(1).toarray())/MEG16.shape[0]
8.0
See also
pyfaust.factparams.ParamsHierarchicalRectMat

4. Simplified Parameters for the Discrete Fourier Transform Factorization

>>> from pyfaust import dft
>>> from pyfaust.fact import hierarchical
>>> from numpy.linalg import norm
>>> # generate a DFT matrix of size 32x32
>>> FDFT = dft(32)
>>> DFT = FDFT.toarray()
>>> # factorize it
>>> FDFT2 = hierarchical(DFT, 'dft');
>>> # test the relative error
>>> (FDFT-FDFT2).norm('fro')/FDFT.norm('fro') # the result is about 1e-16, the factorization is accurate
1.5867087530782683e-16
>>> FDFT

Faust size 32x32, density 0.34375, nnz_sum 352, 6 factor(s):

  • FACTOR 0 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 1 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 2 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 3 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 4 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 5 (complex) SPARSE, size 32x32, density 0.03125, nnz 32
>>> FDFT2

Faust size 32x32, density 0.34375, nnz_sum 352, 6 factor(s):

  • FACTOR 0 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 1 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 2 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 3 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 4 (complex) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 5 (complex) SPARSE, size 32x32, density 0.03125, nnz 32
See also
pyfaust.factparams.ParamsHierarchicalDFT

5. Simplified Parameters for Hadamard Factorization without residual constraints

This factorization parameterization is the same as the one shown in 2. except that there is no constraints at all on residual factors. See pyfaust.factparams.ParamsHierarchicalNoResCons and pyfaust.factparams.ParamsHierarchicalWHTNoResCons for more details.

>>> from pyfaust import wht
>>> from pyfaust.fact import hierarchical
>>> from numpy.linalg import norm
>>> # generate a Hadamard Faust of size 32x32
>>> FH = wht(32)
>>> H = FH.toarray() # the full matrix version
>>> # factorize it
>>> FH2 = hierarchical(H, 'hadamard_simple', backend=2020);
>>> # test the relative error
>>> (FH-FH2).norm('fro')/FH.norm('fro') # the result is about 1e-16, the factorization is accurate
7.85046215906392e-16
>>> FH2

Faust size 32x32, density 0.3125, nnz_sum 320, 5 factor(s):

  • FACTOR 0 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 1 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 2 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 3 (double) SPARSE, size 32x32, density 0.0625, nnz 64
  • FACTOR 4 (double) SPARSE, size 32x32, density 0.0625, nnz 64


See also
pyfaust.factparams.ParamsHierarchicalWHTNoResCons

6. Simplified Parameters for a Rectangular Matrix Factorization (the BSL demo MEG matrix) without residual constraints

The factorization parameterization shown here is the same as in 3. except that there is no constraint at all on residual factors. See pyfaust.factparams.ParamsHierarchicalNoResCons and pyfaust.factparams.ParamsHierarchicalRectMatNoResCons for more details. In the example below the MEG matrix is factorized according to the parameterization shown in 3. (aka "MEG") and on the other hand with the parameterization of interest here (aka "MEG_SIMPLE", with no residual constraints), the approximate accuracy is quite the same so we can conclude that on this case (as in 5.) removing residual constraints can not only simplify the parameterization of hierarchical PALM4MSA but also be as efficient.

>>> from scipy.io import loadmat
>>> from pyfaust import Faust
>>> from pyfaust.fact import hierarchical
>>> from pyfaust.demo import get_data_dirpath
>>> MEG = loadmat(get_data_dirpath()+'/matrix_MEG.mat')['matrix'].T
>>> F1 = hierarchical(MEG, ['MEG', 5, 10, 8], backend=2020)
>>> F2 = hierarchical(MEG, ['MEG_SIMPLE', 5, 10, 8], backend=2020)
>>> # compare the errors:
>>> (F2 - MEG).norm() / Faust(MEG).norm()
0.13033595653237676
>>> (F1 - MEG).norm() / Faust(MEG).norm()
0.12601709513741005


See also
pyfaust.factparams.ParamsHierarchicalRectMatNoResCons
pyfaust.factparams.ParamsHierarchicalRectMat, [1] Le Magoarou L. and Gribonval R., “Flexible multi-layer sparse approximations of matrices and applications”, Journal of Selected Topics in Signal Processing, 2016.

◆ hierarchical_mhtp()

def pyfaust.fact.hierarchical_mhtp (   M,
  hierar_p,
  mhtp_p,
  ret_lambda = False,
  ret_params = False,
  on_gpu = False 
)

Runs the MHTP-PALM4MSA hierarchical factorization algorithm on the matrix M.

This algorithm uses the MHTP-PALM4MSA (pyfaust.fact.palm4msa_mhtp) instead of only PALM4MSA as pyfaust.fact.hierarchical.

Parameters
M(np.ndarray) the numpy array to factorize. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
p(pyfaust.factparams.ParamsHierarchical) is a set of hierarchical factorization parameters. See
hierarchical.
mhtp_p(pypyfaust.factparams.MHTPParams) the instance to define the MHTP algorithm parameters.
on_gpu(bool) if True the GPU implementation is executed.
ret_lambda(bool) set to True to ask the function to return the scale factor (False by default).
ret_params(bool) set to True to ask the function to return the pyfaust.factparams.ParamsHierarchical instance used (False by default).
Returns
See hierarchical

Examples

>>> from pyfaust.fact import hierarchical_mhtp
>>> from pyfaust.factparams import ParamsHierarchical, StoppingCriterion
>>> from pyfaust.factparams import MHTPParams
>>> from pyfaust.proj import sp, normcol, splin
>>> import numpy as np
>>> M = np.random.rand(500, 32)
>>> fact_cons = [splin((500, 32), 5), sp((32,32), 96), sp((32,32), 96)]
>>> res_cons = [normcol((32,32), 1), sp((32,32), 666), sp((32,32), 333)]
>>> stop_crit1 = StoppingCriterion(num_its=200)
>>> stop_crit2 = StoppingCriterion(num_its=200)
>>> # 50 iterations of MHTP will run every 100 iterations of PALM4MSA (each time PALM4MSA is called by the hierarchical algorithm)
>>> mhtp_param = MHTPParams(num_its=50, palm4msa_period=100)
>>> param = ParamsHierarchical(fact_cons, res_cons, stop_crit1, stop_crit2)
>>> param.is_verbose = True
>>> F = hierarchical_mhtp(M, param, mhtp_param)
See also
hierarchical,
Returns
palm4msa_mhtp

◆ palm4msa()

def pyfaust.fact.palm4msa (   M,
  p,
  ret_lambda = False,
  backend = 2016,
  on_gpu = False 
)

Factorizes the matrix M with Palm4MSA algorithm using the parameters set in p.

Parameters
M(np.ndarray) the numpy array to factorize. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
p(pyfaust.factparams.ParamsPalm4MSA) parameters instance to define the algorithm parameters.
ret_lambda(bool) set to True to ask the function to return the scale factor (False by default).
backend(str) the C++ implementation to use (default to 2016, 2020 backend should be faster for most of the factorizations).
on_gpu(bool) if True the GPU implementation is executed (this option applies only to 2020 backend).
Note
If backend==2020 and independently to the StoppingCriterion defined in p, it is possible to stop the algorithm manually at any iteration by the key combination CTRL-C. The last Faust computed in the factorization process will be returned. A typical use case is when the verbose mode is enabled and you see that the error doesn't change anymore or only slightly, you might stop iterations by typing CTRL-C.
Returns
The Faust object resulting of the factorization. if ret_lambda == True then the function returns a tuple (Faust, lambda).

Examples

>>> from pyfaust.fact import palm4msa
>>> from pyfaust.factparams import ParamsPalm4MSA, ConstraintList, StoppingCriterion
>>> import numpy as np
>>> M = np.random.rand(500, 32)
>>> cons = ConstraintList('splin', 5, 500, 32, 'normcol', 1.0, 32, 32)
>>> # or alternatively using pyfaust.proj
>>> # from pyfaust.proj import splin, normcol
>>> # cons = [ splin((500,32), 5), normcol((32,32), 1.0)]
>>> stop_crit = StoppingCriterion(num_its=200)
>>> param = ParamsPalm4MSA(cons, stop_crit)
>>> F = palm4msa(M, param)
>>> F

Faust size 500x32, density 0.22025, nnz_sum 3524, 2 factor(s):

  • FACTOR 0 (double) SPARSE, size 500x32, density 0.15625, nnz 2500
  • FACTOR 1 (double) DENSE, size 32x32, density 1, nnz 1024
See also
pyfaust.factparams.ParamsPalm4msaWHT to factorize a Hadamard matrix using the pyfaust.proj.skperm projector.

◆ palm4msa_mhtp()

def pyfaust.fact.palm4msa_mhtp (   M,
  palm4msa_p,
  mhtp_p,
  ret_lambda = False,
  on_gpu = False 
)

Runs the MHTP-PALM4MSA algorithm to factorize the matrix M.

MHTP stands for Multilinear Hard Tresholding Pursuit. This is a generalization of the Bilinear HTP algorithm describe in [1].

[1] Quoc-Tung Le, Rémi Gribonval. Structured Support Exploration For Multilayer Sparse Matrix Factorization. ICASSP 2021 - IEEE International Conference on Acoustics, Speech and Signal Processing, Jun 2021, Toronto, Ontario, Canada. pp.1-5. hal-03132013

Parameters
M(np.ndarray) the numpy array to factorize. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
palm4msa_p(pyfaust.factparams.ParamsPalm4MSA) instance of parameters to define the PALM4MSA algorithm parameters.
mhtp_p(pyfaust.factparams.MHTPParams) instance of parameters to define the MHTP algorithm parameters.
ret_lambda(bool) set to True to ask the function to return the scale factor (False by default).
on_gpu(bool) if True the GPU implementation is executed.
Returns
  • The Faust object resulting of the factorization.
  • if ret_lambda == True then the function returns a tuple (Faust, lambda).


Example
>>> from pyfaust.fact import palm4msa_mhtp
>>> from pyfaust.factparams import ParamsPalm4MSA, StoppingCriterion, MHTPParams
>>> import numpy as np
>>> from pyfaust.proj import splin, normcol
>>> np.random.seed(42) # just for reproducibility
>>> M = np.random.rand(500, 32)
>>> cons = [ splin((500,32), 5), normcol((32,32), 1.0)]
>>> stop_crit = StoppingCriterion(num_its=200)
>>> param = ParamsPalm4MSA(cons, stop_crit)
>>> # MHTP will run every 100 iterations of PALM4MSA (that is 2 times) for 50 iterations on each factor
>>> mhtp_param = MHTPParams(num_its=50, palm4msa_period=100)
>>> G = palm4msa_mhtp(M, param, mhtp_param)
>>> G
Faust size 500x32, density 0.22025, nnz_sum 3524, 2 factor(s):
  • FACTOR 0 (double) SPARSE, size 500x32, density 0.15625, nnz 2500
  • FACTOR 1 (double) DENSE, size 32x32, density 1, nnz 1024
See also
pyfaust.factparams.MHTPParams, palm4msa pyfaust.fact.hierarchical_mhtp

◆ pinvtj()

def pyfaust.fact.pinvtj (   M,
  nGivens = None,
  tol = 0,
  err_period = 100,
  relerr = True,
  nGivens_per_fac = None,
  enable_large_Faust = False,
**  kwargs 
)

Computes the pseudoinverse of M using svdtj.

Parameters
M(np.ndarray) the matrix to compute the pseudoinverse.
nGivens(int or tuple(int,int)) see svdtj
tol(float) see svdtj (here the error is computed on U.H S^+ V).
err_period(int) see svdtj.
relerr(bool) see svdtj.
nGivens_per_fac(int or tuple(int,int)) see svdtj.
enable_large_Faust(bool) see svdtj.
Returns
The tuple V,Sp,Uh: such that V*numpy.diag(Sp)*Uh is the approximate of M^+.
  • (np.array vector) Sp the inverses of the min(m, n) nonzero singular values in ascending order. Note however that zeros might occur if M is rank r < min(*M.shape).
  • (Faust objects) V, Uh orthonormal Fausts.


Example
>>> from pyfaust.fact import pinvtj
>>> from scipy.sparse import spdiags
>>> import numpy as np
>>> from numpy.random import rand, seed
>>> seed(42)
>>> M = np.random.rand(128, 64)
>>> V, Sp, Uh = pinvtj(M, tol=1e-3)
>>> scipy_Sp = spdiags(Sp, [0], M.shape[1], M.shape[0])
>>> err = np.linalg.norm(V @ scipy_Sp @ Uh @ M - np.eye(64, 64)) / np.linalg.norm(np.eye(64, 64))
>>> print(np.round(err, decimals=6))
0.00012
See also
svdtj

◆ svdtj()

def pyfaust.fact.svdtj (   M,
  nGivens = None,
  tol = 0,
  err_period = 100,
  relerr = True,
  nGivens_per_fac = None,
  enable_large_Faust = False,
**  kwargs 
)

Performs a singular value decomposition and returns the left and right singular vectors as Faust transforms.

Note
this function is based on fact.eigtj which relies on the truncated Jacobi algorithm, hence the 'tj' in the name. See below the examples for further details on how svdtj is defined using eigtj. It's noteworthy that svdtj is still in experimental status, for example it exists cases for which it won't converge (to an abitrary precision, see tol argument). Another important thing is that it needs to compute M.M^H and M^H.M which is not particularly advisable when M is large and dense.
Parameters
M(np.ndarray or scipy.sparse.csr_matrix) a real or complex matrix . The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
nGivens(int or tuple(int, int)) defines the number of Givens rotations that will be used at most to compute U and V. If it is an integer, it will apply both to U and V. If it is a tuple of two integers as nGivens = (JU, JV), JU will be the limit number of rotations for U and JV the same for V. nGivens argument is optional if tol is set but becomes mandatory otherwise.
tol(float) this is the error target on U S V' against M. if error <= tol, the algorithm stops. See relerr below for the error formula. This argument is optional if nGivens is set, otherwise it becomes mandatory.
err_period(int) it defines the period, in number of factors of U or V, the error is compared to tol (reducing the period spares some factors but increases slightly the computational cost because the error is computed more often).
relerr(bool) defines the type of error used to evaluate the stopping criterion defined by tol. true for a relative error ( \(\| U' S V - M\|_F \over \| M \|_F\)) otherwise this is an absolute error ( \(\| U' S V - M\|_F\)).
nGivens_per_fac(int or tuple(int, int)) this argument is the number of Givens rotations to set at most by factor of U and V. If this is an integer it will be the same number of U and V. Otherwise, if it is a tuple of integers (tU, tV), tU will be the number of Givens rotations per factor for U and tV the same for V. By default, this parameter is maximized for U and V, i.e. tU = M.shape[0] / 2, tV = M.shape[1] / 2.
enable_large_Faustsee eigtj.
Note
In order to speed up the error computation of the approximate Faust, the algorithm follows a two-times strategy:
  1. In the first iterations a rough error but less costly error is computed \(\| M \|^2_F - \| S \|^2_F\)
  2. In the latest iterations the precise error (relative or absolute) indicated above (see relerr argument) is computed to reach the targeted accuracy. Don't forget that the err_period argument determines the frequency at which the error is calculated, you might decrease its value in order to obtain an error closer to what you asked. In last resort you can disable the two-times strategy by setting the environment related variable like this (defaulty the value is '0', which means that the two-times strategy is used).
    os.environ['SVDTJ_ALL_TRUE_ERR'] = '1'
Returns
The tuple U,S,V: such that U*numpy.diag(S)*V.H is the approximate of M.
  • (np.array vector) S the singular values in any order.
  • (Faust objects) U,V orthonormal Fausts.
Examples
>>> from pyfaust.fact import svdtj
>>> from numpy.random import rand, seed
>>> import numpy as np
>>> from scipy.sparse import spdiags
>>> seed(42) # just for reproducibility
>>> M = rand(16, 32)
>>> # Factoring by specifying the number of Givens rotations
>>> U1, S1, V1 = svdtj(M, 4096, enable_large_Faust=True)
>>> S1_ = spdiags(S1, [0], U1.shape[0], V1.shape[0])
>>> np.allclose(U1 @ S1_ @ V1.H, M)
True
>>> # Specifying a different number of rotations for U and V
>>> # Because U is smaller it should need less rotations
>>> U2, S2, V2 = svdtj(M, (2400, 3200), enable_large_Faust=True)
>>> S2_ = spdiags(S2, [0], U2.shape[0], V2.shape[0])
>>> np.allclose(U2 @ S2_ @ V2.H, M)
True
>>> # Factoring according to an approximate accuracy target
>>> U3, S3, V3 = svdtj(M, tol=1e-12, enable_large_Faust=False)
>>> S3_ = spdiags(S3, [0], U3.shape[0], V3.shape[0])
>>> np.allclose(U3 @ S3_ @ V3.H, M)
True
>>> # verify the relative error is lower than 1e-12
>>> np.linalg.norm(U3 @ S3_ @ V3.H - M) / np.linalg.norm(M)
9.122443356294125e-13
>>> # try with an absolute tolerance (the previous one was relative to M)
>>> U4, S4, V4 = svdtj(M, tol=1e-6, relerr=False, enable_large_Faust=False)
>>> S4_ = spdiags(S4, [0], U4.shape[0], V4.shape[0])
>>> np.allclose(U4 @ S4_ @ V4.H, M)
True
>>> # verify the absolute error is lower than 1e-6
>>> np.linalg.norm(U4 @ S4_ @ V4.H - M)
1.2044207331117425e-11
>>> # try a less accurate approximate to get less factors
>>> U5, S5, V5 = svdtj(M, nGivens=(256, 512), tol=1e-1, relerr=True, enable_large_Faust=False)
>>> S5_ = spdiags(S5, [0], U5.shape[0], V5.shape[0])
>>> # verify the absolute error is lower than 1e-1
>>> np.linalg.norm(U5 @ S5_ @ V5.H - M) / np.linalg.norm(M)
0.09351811486725303
>>>
>>> len(V1) # it should be 4096 / nGivens_per_fac, which is (M.shape[1] // 2) = 256
256
>>> len(U1) # it should be 4096 / nGivens_per_fac, which is (M.shape[0] // 2) = 512
100
>>> # but it is not, svdtj stopped automatically extending U1 because the error stopped enhancing
>>> # (it can be verified with verbosity=1)
>>> (len(U2), len(V2))
(100, 200)
>>> (len(U3), len(V3))
(64, 256)
>>> (len(U4), len(V4))
(64, 200)
>>> # not surprisingly U5 and V5 use the smallest number of factors (nGivens and tol were the smallest)
>>> (len(U5), len(V5))
(32, 32)
Explanations

If we call svdtj on the matrix M, it makes two internal calls to eigtj. In Python it would be:

  1. D1, W1 = eigtj(M.dot(M.H), next_args...)
  2. D2, W2 = eigtj(M.H.dot(M), next_args...)

It gives the following equalities (ignoring the fact that eigtj computes approximations):

\[W_1 D_1 W_1^* = M M^*\]

\[W_2 D_2 W_2^* = M^* M\]

But because of the SVD \( M = USV^* \) we also have:

\[MM^* = U S V^* V S U^* = U S^2 U^* = W_1 D_1 W_1^*\]

\[M^* M = V S U^* U S V^* = V S^2 V^* = W_2 D_2 W_2^*\]

It allows to identify the left singular vectors of M to W1, and likewise the right singular vectors to W2.

To compute a consistent approximation of S we observe that U and V are orthogonal hence \(S = U^* M V\) so we ignore the off-diagonal coefficients of the approximation and take \(S = diag(U^* M V) \approx diag(W_1^* M W_2)\)

See also
eigtj
pyfaust.demo
The pyfaust demo module.
Definition: demo.py:1
pyfaust.factparams.ParamsHierarchicalRectMatNoResCons
This parameter class is the same as ParamsHierarchicalRectMat except that there is no residual factor...
Definition: factparams.py:2049
pyfaust.fact.pinvtj
def pinvtj(M, nGivens=None, tol=0, err_period=100, relerr=True, nGivens_per_fac=None, enable_large_Faust=False, **kwargs)
Computes the pseudoinverse of M using svdtj.
Definition: fact.py:345
pyfaust.norm
def norm(F, ord='fro', **kwargs)
Returns Faust.norm(F, ord)` ornumpy.linalg.norm(F, ord)`` depending of F type.
Definition: __init__.py:3917
pyfaust.seed
def seed(s)
(Re)Initializes the pyfaust pseudo-random generator.
Definition: __init__.py:5334
pyfaust.dft
def dft(n, normed=True, dev='cpu', diag_opt=False)
Constructs a Faust F implementing the Discrete Fourier Transform (DFT) of order n.
Definition: __init__.py:4293
pyfaust.rand
def rand(num_rows, num_cols, num_factors=None, dim_sizes=None, density=None, fac_type='sparse', per_row=True, dev='cpu', dtype='float64', field=None, seed=0)
Generates a random Faust.
Definition: __init__.py:5016
pyfaust.fact.svdtj
def svdtj(M, nGivens=None, tol=0, err_period=100, relerr=True, nGivens_per_fac=None, enable_large_Faust=False, **kwargs)
Performs a singular value decomposition and returns the left and right singular vectors as Faust tran...
Definition: fact.py:227
pyfaust.fact.palm4msa
def palm4msa(M, p, ret_lambda=False, backend=2016, on_gpu=False)
Factorizes the matrix M with Palm4MSA algorithm using the parameters set in p.
Definition: fact.py:567
pyfaust.fact.hierarchical_mhtp
def hierarchical_mhtp(M, hierar_p, mhtp_p, ret_lambda=False, ret_params=False, on_gpu=False)
Runs the MHTP-PALM4MSA hierarchical factorization algorithm on the matrix M.
Definition: fact.py:711
pyfaust.factparams.ParamsHierarchicalSquareMat
Definition: factparams.py:1885
pyfaust.fact.eigtj
def eigtj(M, nGivens=None, tol=0, err_period=100, order='ascend', relerr=True, nGivens_per_fac=None, verbosity=0, enable_large_Faust=False)
Performs an approximate eigendecomposition of M and returns the eigenvalues in W along with the corre...
Definition: fact.py:470
pyfaust.factparams.ParamsHierarchicalDFT
The simplified parameterization class for factorizing a DFT matrix using the hierarchical factorizati...
Definition: factparams.py:1689
pyfaust.proj
This module provides matrix projectors.
Definition: proj.py:1
pyfaust.fact.palm4msa_mhtp
def palm4msa_mhtp(M, palm4msa_p, mhtp_p, ret_lambda=False, on_gpu=False)
Runs the MHTP-PALM4MSA algorithm to factorize the matrix M.
Definition: fact.py:658
pyfaust.wht
def wht(n, normed=True, dev="cpu", dtype='float64')
Constructs a Faust implementing the Walsh-Hadamard Transform (WHT) of order n.
Definition: __init__.py:4142
pyfaust.factparams
The module for the parameterization of FAuST's algorithms (Palm4MSA and Hierarchical Factorization)....
Definition: factparams.py:1
pyfaust.factparams.ParamsHierarchicalWHTNoResCons
The simplified parameterization class for factorizing a Hadamard matrix with the hierarchical factori...
Definition: factparams.py:1841
pyfaust.fact.hierarchical
def hierarchical(M, p, ret_lambda=False, ret_params=False, backend=2016, on_gpu=False)
Factorizes the matrix M using the hierarchical PALM4MSA algorithm according to the parameters set in ...
Definition: fact.py:1014
pyfaust.demo.get_data_dirpath
def get_data_dirpath(silent=True)
Returns the data directory path which varies according to the way pyfaust was installed.
Definition: demo.py:60
pyfaust.fact.butterfly
def butterfly(M, type="bbtree", perm=None, diag_opt=False, mul_perm=None)
Factorizes M according to a butterfly support and optionally a permutation using the algorithms descr...
Definition: fact.py:1245
pyfaust.fact
The pyfaust factorization module.
Definition: fact.py:1