[PYTHON] Sparse modeling: 2.3 Construction of Glasman matrix

I tried to build the Glasman matrix in Chapter 2.3 of Sparse Modeling by Michael Elad and translated by Toru Tamaki in Python.

In sparse modeling $(P_0): \min_x ||\mathbf{x}||_0 \, \rm{subject \, to} \, \mathbf{b}=\mathbf{Ax}$ To solve the problem. Let $ \ mathbf {A} $ be the $ n × m $ matrix $ (m> n) $.

Here, the closer $ \ mathbf {A} $ is to the orthogonal matrix (the closer $ \ mathbf {A} ^ T \ mathbf {A} $ is to $ \ mathbf {I} $), the closer to obtain the solution of $ P_0 $. Will be easier.

Since $ \ mathbf {A} $ is not a square matrix, the off-diagonal elements of $ \ mathbf {A} ^ T \ mathbf {A} $ cannot all be 0, and its lower limit is

\mu(\mathbf{A}) \geq \sqrt{\frac{m-n}{n(m-1)}}

It will be. This $ \ mu (\ mathbf {A}) $ is called __ mutual coherence __ of the matrix $ \ mathbf {A} $.

When $ \ mathbf {x} $ is found as a solution candidate for $ P_0 $ by some algorithm, $ \ mathbf {x} $ is

||\mathbf{x}||_0 < \frac{1}{2}(1+\frac{1}{\mu(\mathbf{A})})

It can be said that $ \ mathbf {x} $ is the only sparse solution under the constraint of $ \ mathbf {b} = \ mathbf {Ax} $ if the condition of is satisfied (Theorem 2.5). .. So if you're looking for a sparse $ \ mathbf {x} $ that satisfies $ \ mathbf {b} = \ mathbf {Ax} $ with some approximate algorithm, the value of $ \ mu (\ mathbf {A}) $ It can be said that the smaller is, the easier it is to find the answer on a looser basis.

That's why the lower the mutual coherence of $ \ mathbf {A} $, the better.

\mu(\mathbf{A}) = \sqrt{\frac{m-n}{n(m-1)}}

A matrix whose mutual coherence matches the lower limit is called the Glasman matrix. It can be said that the Glasman matrix is a $ n × m $ matrix that is as close to an orthogonal matrix as possible.

I think that it would be useful to solve the $ P_0 $ problem if such a Glasman matrix can be constructed, but the Glasman matrix is calculated numerically in "2.3 Construction of the Glasman matrix" of sparse modeling. I had the MATLAB code to do, so I ported it to Python (with a long introduction).

The algorithm itself shrinks a large value of the Gramian matrix $ \ mathbf {G} = \ mathbf {A} ^ T \ mathbf {A} $ and then performs SVD conversion so that the rank of the matrix is N. When iterated, a Gram matrix $ \ mathbf {G} $ whose off-diagonal component value was close to the lower limit was obtained, and finally, the Gram matrix $ \ mathbf {G} $ which was iteratively processed by SVD conversion became a Glasman matrix. You get $ \ hat {\ mathbf {A}} $.

See ipynb on Github below for Python code to build the Glasman matrix. https://github.com/kibo35/sparse-modeling/blob/master/ch02.ipynb

It is suspicious to get $ \ hat {\ mathbf {A}} $ which became the Glasman matrix at the end, but I think that the processing to the Gram matrix gives the same result as the textbook.

The upper left is the gram matrix $ \ mathbf {G} = \ mathbf {A} ^ T \ mathbf {A} $ of the redundant system matrix $ \ mathbf {A} $ (50 rows x 100 columns) given by the normal distribution. I will. The upper right is the Gram matrix after 10000 iterations.

If you look at the image with the absolute value in the lower row, you can see that the values of the off-diagonal components are uniform.

The figure above plots the values of the off-diagonal components of the Gram matrix at each iteration.

Anyway, from Chapter 2 of sparse modeling, we can see that the mutual coherence of the redundant system matrix $ \ mathbf {A} $ $ \ mu (\ mathbf {A}) $ (maximum absolute value of the off-diagonal components of the Gram matrix). If you know the value of) and find a solution $ \ mathbf {x} $ whose norm is smaller than a certain threshold value, you can say that it is a unique solution, so do your best. Let's solve the $ P_0 $ problem.