This notebook put the focus on the proj module of the pyfaust API. This module provides a bunch of projectors which are for instance necessary to implement proximal operators in PALM algorithms.
Indeed these projectors matter in the parameterization of the PALM4MSA and hierarchical factorization algorithms, so let's maintain their configuration as simple as possible by using projectors!
First let's explain some generalities about projectors:
proj_gen.The general pattern to use a projector unfolds in two steps:
Let's see how to define and use the projectors in the code. For the brief math definitions, I let you consult this document. Remember that the projector API is documented too, you'll find the link for each projector below. Last, if you're looking for a reference about proximal operators here it is: proximity operators.
This projector performs a projection onto matrices with a prescribed sparsity. It governs the global sparsity of a matrix given an integer k. The matrix $ A \in \mathbb{R}^{m \times n }, A = (a_{ij}), {0<= i < m, 0<= j < n}$, is projected to the closest matrix $ B = (b_{ij})$ such that $ \|B\|_0 = \#\{(i,j): b_{ij} \neq 0 \} \leq k$ which implies , if $k < mn$, that some entries of $A$ are kept in $B$ and others are set to zero. The projector keeps the $k$ most significant values (in term of absolute values or magnitude).
Let's try on an example, here a random matrix.
from numpy.random import rand
A = rand(5,5)*100
print("A=\n", A)
A= [[44.09291282 69.65449995 79.01609433 21.15468406 29.01133231] [36.17382938 4.5085547 82.54968789 2.97148275 28.46251655] [69.25452372 68.16201213 85.65208859 15.03199386 71.86439886] [84.70875566 18.00189631 73.66586237 93.34860797 80.63185061] [78.17918182 35.39168531 97.00055955 7.8148175 37.27405543]]
import pyfaust
from pyfaust.proj import sp
# 1. instantiate the projector
k = 2
p = sp(A.shape, k, normalized=False)
# 2. project the matrix through it
B = p(A)
print("B=\n", B)
B= [[ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 93.34860797 0. ] [ 0. 0. 97.00055955 0. 0. ]]
The projector is simply defined by the input matrix shape and the integer k to specify the targeted sparsity.
Optional normalization:
As you noticed, the argument normalized is set to False in the projector definition. This is the default behaviour. When normalized is True, the result $B$ is normalized according to its Frobenius norm.
The next example gives you a concrete view of what happens when normalized is True.
from numpy.linalg import norm
from numpy import allclose
pnorm = sp(A.shape, k, normalized=True)
C = pnorm(A)
print("B/norm(B, 'fro') == C: ", allclose(B/norm(B,'fro'),C))
B/norm(B, 'fro') == C: True
Sparsity and optional positivity:
It is also possible to "filter" the negative entries of A by setting the pos argument of sp to True.
You can see the projector as a pipeline, the first stage is to filter out the negative values, then the sp projector is applied and finally the resulting image is normalized if normalized==True.
The following example shows how the projector operates depending on combinations of pos and normalized values.
p_pos = sp(A.shape, k, normalized=False, pos=True)
print(p_pos(A))
[[ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 93.34860797 0. ] [ 0. 0. 97.00055955 0. 0. ]]
Well, it's exactly the same as the In [2] output. The reason is quite obvious, it's because A doesn't contain any negative value, so let's try on a copy of A where we set the p_pos(A) nonzeros to negative values.
from numpy import nonzero
D = A.copy()
D[nonzero(p_pos(A))] *= -1
print(p_pos(D))
[[ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 85.65208859 0. 0. ] [84.70875566 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ]]
The entries selected when p_pos(A) is applied are now skipped because they are negative in D, so p_pos(D) selects the new greatest two values of A in term of magnitude.
What happens now if all values of the matrix are negative ? Let's see it in the next example.
E = - A.copy()
print(p_pos(E))
[[0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.]]
That should not be surprising that the resulting matrix is a zero matrix, indeed E contains only negative values which are all filtered by setting pos=True in the p_pos definition.
A last question remains: what would happen if we normalize the output matrix when pos==True and the input matrix is full of negative values ?
The response is simple: a division by zero error would be raised because the norm of a zero matrix is zero, hence it's not possible to normalize.
They are very similar to the sp projector except that splin governs the integer sparsity on a row basis and spcol does it by columns as indicated by the suffix name.
Look at the two short examples, just to be sure.
from pyfaust.proj import splin, spcol
pl = splin(A.shape, k) # reminder: k == 2
pc = spcol(A.shape, k)
B1 = pl(A)
B2 = pc(A)
print("B1=\n", B1)
print("\bB2=\n", B2)
B1= [[ 0. 69.65449995 79.01609433 0. 0. ] [36.17382938 0. 82.54968789 0. 0. ] [ 0. 0. 85.65208859 0. 71.86439886] [84.70875566 0. 0. 93.34860797 0. ] [78.17918182 0. 97.00055955 0. 0. ]] B2= [[ 0. 69.65449995 0. 21.15468406 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 68.16201213 85.65208859 0. 71.86439886] [84.70875566 0. 0. 93.34860797 80.63185061] [78.17918182 0. 97.00055955 0. 0. ]]
Here the k most significant values are chosen (by rows for splin or by columns for spcol) and the image normalization is disabled.
As for the SP projector, it is possible to incorporate a normalization and/or positivity constraint passing normalized=True and pos=True to the functor constructor.
The projector splincol tries to constrain the sparsity both by columns and by rows and I wrote it "tries" because there is not always a solution. The use is again the same.
from pyfaust.proj import splincol
plc = splincol(A.shape, k)
print(plc(A))
[[ 0. 69.65449995 79.01609433 21.15468406 0. ] [36.17382938 0. 82.54968789 0. 0. ] [ 0. 68.16201213 85.65208859 0. 71.86439886] [84.70875566 0. 0. 93.34860797 80.63185061] [78.17918182 0. 97.00055955 0. 0. ]]
The image matrix support is in fact the union set of the supports obtained through splin and spcol projectors (that's the reason why there is not always a solution). You can refer to this documentation page which demonstrates in an example how this union is defined.
Another projector for the same purpose but more precise is available for square matrices only. It is named skperm, you'll find its API doc here. In brief, it is based on a derivate of the hungarian algorithm.
Another matrix projector is the blockdiag projector. As its name suggests it projects onto the closest block-diagonal matrix with a prescribed structure.
The block-diagonal structure can be defined by the list of the shapes of the block diagonal submatrices you want to keep from the input matrix into the output matrix.
An example will show how it works easily:
%matplotlib inline
from pyfaust.proj import blockdiag
from pyfaust import Faust
R = rand(15,25)
pbd = blockdiag(R.shape, [(1,1), (1,12) ,(R.shape[0]-2, R.shape[1]-13)]) # shapes of blocks in second argument
# show it as a Faust composed of a single factor
Faust([pbd(R)]).imshow()
This blockdiag projector above is defined in order to keep three blocks of the input matrix A, from the upper-left to the lower-right: the first block is the singleton block composed only of the entry (0,0), the second block is a bit of the next row starting from entry (1,1) and finishing to the entry (1, 12) (its shape is (1,12)) and the final block starts from the element (2,13) to finish on the element (R.shape[0]-1, R.shape[1]-1). It's important that the list of blocks covers the whole matrix from its entry (0,0) to its entry (R.shape[0]-1, R.shape[1]-1) or the projector will end up on an error. In other words, if you sum together the first (resp. the second) coordinates of all pair of shapes you must find R.shape[0] (resp. R.shape[1]).
These projectors all return the closest corresponding structured matrix (circ being the short name for circulant). For detailed definitions of these kinds of matrices you can refer to wikipedia:
The output is constant along each 'diagonal' (resp. 'anti-diagonal'), where the corresponding constant value is the mean of the values of the input matrix along the same diagonal (resp. anti-diagonal).
In the following example, a Faust is constructed with in first factor a circulant matrix, in second position a toeplitz matrix and at the end a hankel matrix.
from pyfaust.proj import circ, toeplitz, hankel
CI = rand(10,10) # circ proj input
TI = rand(10,15) # toeplitz proj input
HI = rand(15,10) # hankel proj input
cp = circ(CI.shape)
tp = toeplitz(TI.shape)
hp = hankel(HI.shape)
F = Faust([cp(CI), tp(TI), hp(HI)])
F.imshow()
You should recognize clearly the structures of a circulant matrix, a toeplitz matrix and a hankel matrix.
Note that these projectors are also capable to receive the normalized and pos keyword arguments we've seen before.
The API documentation will give you other examples:
The pyfaust.proj module provides two projectors, normlin (resp. normcol) that projects a matrix onto the closest matrix with rows (resp. columns) of a prescribed 2-norm.
Let's try them:
from pyfaust.proj import normcol, normlin
from numpy.linalg import norm
pnl = normlin(A.shape, .2)
pnc = normcol(A.shape, .2)
# let's verify the norm for one column obtained by normlin
print(norm(pnl(A)[2,:]))
# and the same for normcol
print(norm(pnc(A)[:,2]))
0.2 0.2
Something that is important to notice is the particular case of zero columns or rows. When the NORMLIN (resp. NORMCOL) projector encounters a zero row (resp. a zero column) it simply ignores it. Let's try:
B = A.copy()
B[:,2] = 0
print(pnc(B)[:,2])
[0. 0. 0. 0. 0.]
The column 2 is set to zero in B and stays to zero in the projector image matrix.
The supp projector projects a matrix onto the closest matrix with a prescribed support. In other words, it preserves the matrix entries lying on this support. The others are set to zero.
The support must be defined by a binary matrix (the dtype must be the same as the input matrix though).
from pyfaust.proj import supp
from numpy import eye
# keep only the diagonal of A
ps = supp(eye(*A.shape)) # by default normalized=False and pos=False
print(ps(A))
[[44.09291282 0. 0. 0. 0. ] [ 0. 4.5085547 0. 0. 0. ] [ 0. 0. 85.65208859 0. 0. ] [ 0. 0. 0. 93.34860797 0. ] [ 0. 0. 0. 0. 37.27405543]]
This projector is also capable to receive the normalized and pos keyword arguments.
This last projector is really simple, it returns a constant matrix whatever is the input matrix. The way it's instantiated is very similar to the SUPP projector. Look at its documentation to get an example: const(ant) projector.
Thanks for reading this notebook, you'll find others on the FAµST website. Any feedback is welcome, all contact information is as well available on the website.
Note: this notebook was executed using the following pyfaust version:
import pyfaust
pyfaust.version()
'3.0.11'