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:

- They are all functor objects (objects that you can call as a function).
- They are all types defined by child classes of the parent abstract class
`proj_gen`

.

The general pattern to use a projector unfolds in two steps:

- Instantiate the projector passing the proper arguments.
- Call this projector (again, as a function) on the matrix you're working on. This step is optional or for test purposes, because generally it is the algorithm implementation responsibility to call the projectors. You just need to feed the algorithms (PALM4MSA for example) with them.

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.

In [1]:

```
from numpy.random import rand
A = rand(5,5)*100
print("A=\n", A)
```

In [2]:

```
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)
```

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.

In [3]:

```
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.

In [4]:

```
p_pos = sp(A.shape, k, normalized=False, pos=True)
print(p_pos(A))
```

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.

In [5]:

```
from numpy import nonzero
D = A.copy()
D[nonzero(p_pos(A))] *= -1
print(p_pos(D))
```

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.

In [6]:

```
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.

In [7]:

```
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)
```

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.

In [8]:

```
from pyfaust.proj import splincol
plc = splincol(A.shape, k)
print(plc(A))
```

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:

In [9]:

```
%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.

In [10]:

```
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:

In [11]:

```
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:

In [12]:

```
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).

In [13]:

```
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))
```

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:

In [14]:

```
import pyfaust
pyfaust.version()
```

Out[14]:

'3.0.11'