# SPGL1 demo

This notebook contains a Python implementation of the original examples from SPGL1 MATLAB solver



In [1]:

%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt

from scipy.sparse import spdiags
from scipy.sparse.linalg import lsqr as splsqr
from spgl1.lsqr import lsqr
from spgl1 import spgl1, spg_lasso, spg_bp, spg_bpdn, spg_mmv
from spgl1.spgl1 import norm_l1nn_primal, norm_l1nn_dual, norm_l1nn_project
from spgl1.spgl1 import norm_l12nn_primal, norm_l12nn_dual, norm_l12nn_project

# Initialize random number generators
np.random.seed(43273289)



## Lasso



In [2]:

# Create random m-by-n encoding matrix and sparse vector
m = 50
n = 128
k = 14
[A,Rtmp] = np.linalg.qr(np.random.randn(n,m),'reduced')
A  = A.T
p  = np.random.permutation(n)
p = p[0:k]
x0 = np.zeros(n)
x0[p] = np.random.randn(k)



Solve the underdetermined LASSO problem for $||x||_1 <= \pi$:

$$min.||Ax-b||_2 \quad subject \quad to \quad ||x||_1 <= \pi$$


In [3]:

b = A.dot(x0)
tau = np.pi
x,resid,grad,info = spg_lasso(A, b, tau, verbosity=1)

print()
print('%s%s%s' % ('-'*35,' Solution ','-'*35))
print('nonzeros(x) = %i,   ||x||_1 = %12.6e,   ||x||_1 - pi = %13.6e' % \
(np.sum(abs(x)>1e-5), np.linalg.norm(x,1), np.linalg.norm(x,1)-np.pi))
print('%s' % ('-'*80))




================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 3.14e+00
Two-norm of b         : 3.40e+00

Optimality tol        : 1.00e-04
Target one-norm of x  : 3.14e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :      500

EXIT -- Optimal solution found

Products with A     :       8        Total time   (secs) :     0.0
Products with A^H   :       8        Project time (secs) :     0.0
Newton iterations   :       0        Mat-vec time (secs) :     0.0
Line search its     :       0        Subspace iterations :       0

----------------------------------- Solution -----------------------------------
nonzeros(x) = 7,   ||x||_1 = 3.141593e+00,   ||x||_1 - pi =  8.881784e-16
--------------------------------------------------------------------------------



## BP

Solve the basis pursuit (BP) problem:

$$min. ||x||_1 \quad subject \quad to \quad Ax = b$$


In [4]:

b = A.dot(x0) # signal

plt.figure()
plt.plot(x,'b')
plt.plot(x0,'ro')
plt.legend(('Recovered coefficients','Original coefficients'))
plt.title('Basis Pursuit');




================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 3.40e+00

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :      500

iterr      Objective   Relative Gap  Rel Error      gnorm   stepg   nnz_x   nnz_g     tau

0  3.3978494e+00  0.0000000e+00   1.00e+00  8.995e-01     0.0       0       0  1.2834830e+01
1  1.4610032e+00  1.6532743e+00   1.00e+00  3.459e-01    -0.6      67       1
2  1.2957367e+00  1.4480527e+00   1.00e+00  3.057e-01     0.0      53       1
3  1.0445582e+00  1.0186987e+00   1.00e+00  2.394e-01     0.0      33       1
4  9.8320803e-01  1.9671347e+00   9.83e-01  3.006e-01     0.0      17       1
5  1.0408537e+00  2.7237524e+00   1.00e+00  3.273e-01     0.0      15       1
6  8.3766484e-01  4.4627639e-01   8.38e-01  1.829e-01     0.0      20       1
7  8.2370803e-01  2.1784189e-01   8.24e-01  1.661e-01     0.0      16       1
8  8.1471442e-01  1.0576174e-01   8.15e-01  1.582e-01     0.0      15       1
9  8.1352564e-01  1.2833205e-01   8.14e-01  1.618e-01     0.0      12       2
10  8.1346284e-01  1.2136051e-01   8.13e-01  1.595e-01     0.0      12       3
12  8.1283336e-01  1.6431186e-03   8.13e-01  1.509e-01     0.0      12      12  1.7214207e+01
20  1.0262800e-01  3.7124664e-02   1.03e-01  1.684e-02     0.0      55       5
30  5.8775070e-02  2.1151854e-02   5.88e-02  9.905e-03     0.0      26       5
39  5.6081993e-02  3.4487799e-03   5.61e-02  9.059e-03     0.0      16      16  1.7561410e+01
40  1.3497461e-02  2.3368157e-02   1.35e-02  2.563e-03     0.0      61      17
50  6.5248570e-03  2.5729401e-02   6.52e-03  2.012e-03    -0.6      44       8
60  3.6779037e-03  8.0294188e-03   3.68e-03  8.458e-04    -0.9      30      30
70  1.7198692e-03  6.8520394e-04   1.72e-03  2.801e-04     0.0      17      17
80  1.3668596e-03  8.5588976e-05   1.37e-03  2.204e-04     0.0      14      14
82  1.3666459e-03  1.8193225e-05   1.37e-03  2.163e-04     0.0      14      14  1.7570044e+01
90  1.3860943e-04  4.2028246e-05   1.39e-04  2.222e-05     0.0      14      14
98  8.6425580e-05  5.4782379e-05   8.64e-05  1.545e-05     0.0      14      14

EXIT -- Found a root

Products with A     :     131        Total time   (secs) :     0.1
Products with A^H   :      99        Project time (secs) :     0.0
Newton iterations   :       4        Mat-vec time (secs) :     0.0
Line search its     :      32        Subspace iterations :       0




In [5]:

plt.figure()
plt.plot(info['xnorm1'], info['rnorm2'], '.-k')
plt.xlabel(r'$||x||_1$')
plt.ylabel(r'$||r||_2$')
plt.title('Sampled Pareto curve')

plt.figure()
plt.plot(np.arange(info['niters']), info['rnorm2']/max(info['rnorm2']), '.-k')
plt.plot(np.arange(info['niters']), info['xnorm1']/max(info['xnorm1']), '.-r')
plt.xlabel(r'#iter')
plt.ylabel(r'$||r||_2 & ||x||_1$');
plt.title('Cost functions');






## BPDN

Solve the basis pursuit denoise (BPDN) problem:

$$min. ||x||_1 \quad subject \quad to \quad ||Ax - b||_2 <= 0.1$$


In [6]:

b = A.dot(x0) + np.random.randn(m) * 0.075
sigma = 0.10  #     % Desired ||Ax - b||_2
x,resid,grad,info = spg_bpdn(A, b, sigma, iter_lim=10, verbosity=2)

plt.figure()
plt.plot(x,'b')
plt.plot(x0,'ro')
plt.legend(('Recovered coefficients','Original coefficients'))
plt.title('Basis Pursuit Denoise');




================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 3.49e+00

Optimality tol        : 1.00e-04
Target objective      : 1.00e-01

Basis pursuit tol     : 1.00e-06
Maximum iterations    :       10

iterr      Objective   Relative Gap  Rel Error      gnorm   stepg   nnz_x   nnz_g     tau

0  3.4876525e+00  0.0000000e+00   9.71e-01  8.998e-01     0.0       0       0  1.3130733e+01
1  1.5278506e+00  1.5445083e+00   9.35e-01  3.512e-01    -0.6      70       1
2  1.3785367e+00  1.3210423e+00   9.27e-01  3.004e-01     0.0      48       1
3  1.1435885e+00  1.0754951e+00   9.13e-01  2.436e-01     0.0      36       1
4  1.0650070e+00  1.8808347e+00   9.06e-01  3.058e-01     0.0      24       1
5  1.1055544e+00  1.9676460e+00   9.10e-01  2.609e-01     0.0      14       1
6  9.7104698e-01  5.0757852e-01   8.71e-01  1.910e-01     0.0      21       1
7  9.5573059e-01  2.5393233e-01   8.56e-01  1.677e-01     0.0      20       1
8  9.5003110e-01  1.6053048e-01   8.50e-01  1.613e-01     0.0      19       1
9  9.4642238e-01  6.3958437e-02   8.46e-01  1.545e-01     0.0      16       1
10  9.4667260e-01  1.0342969e-01   8.47e-01  1.576e-01     0.0      15       1

ERROR EXIT -- Too many iterations

Products with A     :      13        Total time   (secs) :     0.0
Products with A^H   :      11        Project time (secs) :     0.0
Newton iterations   :       1        Mat-vec time (secs) :     0.0
Line search its     :       2        Subspace iterations :       0



## BPDN with non-negative solution

We repeat the same procedure but we have only positive elements in the x vector. We compare spgl1 with L1 norms and with L1NN norms



In [7]:

x0 = np.zeros(n)
x0[p] = np.abs(np.random.randn(k))
b = A.dot(x0) # signal

x,resid,grad,info = spg_bp(A, b, iter_lim=20, verbosity=1)

xnn,residnn,gradnn,infonn = spg_bp(A, b, iter_lim=20, verbosity=1,
project=norm_l1nn_project,
primal_norm=norm_l1nn_primal,
dual_norm=norm_l1nn_dual)

plt.figure()
plt.plot(x,'b')
plt.plot(xnn,'--g')
plt.plot(x0,'ro')
plt.legend(('Recovered coefficients', 'Recovered coefficients NNnorms','Original coefficients'))
plt.title('Basis Pursuit');




================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 1.45e+00

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :       20

ERROR EXIT -- Too many iterations

Products with A     :      24        Total time   (secs) :     0.0
Products with A^H   :      21        Project time (secs) :     0.0
Newton iterations   :       2        Mat-vec time (secs) :     0.0
Line search its     :       3        Subspace iterations :       0

================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 1.45e+00

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :       20

ERROR EXIT -- Too many iterations

Products with A     :      24        Total time   (secs) :     0.0
Products with A^H   :      21        Project time (secs) :     0.0
Newton iterations   :       2        Mat-vec time (secs) :     0.0
Line search its     :       3        Subspace iterations :       0



## BP with complex numbers

Solve the basis pursuit (BP) problem in COMPLEX variables:

$$min. ||z||_1 \quad subject \quad to \quad Az = b$$


In [8]:

from scipy.sparse.linalg import LinearOperator

class partialFourier(LinearOperator):
def __init__(self, idx, n):
self.idx = idx
self.n = n
self.shape = (len(idx), n)
self.dtype = np.complex128
def _matvec(self, x):
# % y = P(idx) * FFT(x)
z = np.fft.fft(x) / np.sqrt(n)
return z[idx]
def _rmatvec(self, x):
z = np.zeros(n,dtype=complex)
z[idx] = x
return np.fft.ifft(z) * np.sqrt(n)

# % Create partial Fourier operator with rows idx
idx = np.random.permutation(n)
idx = idx[0:m]
opA = partialFourier(idx, n)

# % Create sparse coefficients and b = 'A' * z0;
z0 = np.zeros(n,dtype=complex)
z0[p] = np.random.randn(k) + 1j * np.random.randn(k)
b = opA.matvec(z0)

plt.figure()
plt.plot(z.real,'b+',markersize=15.0)
plt.plot(z0.real,'bo')
plt.plot(z.imag,'r+',markersize=15.0)
plt.plot(z0.imag,'ro')
plt.legend(('Recovered (real)', 'Original (real)', 'Recovered (imag)', 'Original (imag)'))
plt.title('Complex Basis Pursuit');




================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 2.78e+00

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :      500

iterr      Objective   Relative Gap  Rel Error      gnorm   stepg   nnz_x   nnz_g     tau

0  2.7830815e+00  0.0000000e+00   1.00e+00  7.891e-01     0.0       0       0  9.8155346e+00
1  1.4897141e+00  1.7264979e+00   1.00e+00  3.852e-01    -0.6      91       0
2  1.3104225e+00  1.1222827e+00   1.00e+00  2.978e-01     0.0      62       0
3  1.1126604e+00  6.4597283e-01   1.00e+00  2.023e-01     0.0      28       0
4  1.0695040e+00  8.0305624e-01   1.00e+00  2.541e-01     0.0      27       0
5  1.0417686e+00  5.7568961e-01   1.00e+00  1.964e-01     0.0      18       0
6  1.0237323e+00  7.8159888e-02   1.00e+00  1.646e-01     0.0      17       0
7  1.0222962e+00  2.9842295e-02   1.00e+00  1.595e-01     0.0      17       0
8  1.0216646e+00  1.9735523e-02   1.00e+00  1.589e-01     0.0      16       0
9  1.0215637e+00  3.4485995e-02   1.00e+00  1.601e-01     0.0      14       0
10  1.0216029e+00  2.7849988e-02   1.00e+00  1.601e-01     0.0      15       0  1.6336186e+01
20  1.4943828e-01  2.2501214e-01   1.49e-01  2.871e-02     0.0      53       0
30  1.0102794e-01  8.8737221e-02   1.01e-01  1.694e-02     0.0      35       2
40  9.3739382e-02  1.2825688e-02   9.37e-02  1.319e-02     0.0      21       2
43  9.3680336e-02  8.4527745e-04   9.37e-02  1.242e-02     0.0      21       6  1.7042723e+01
50  1.1642922e-02  2.3009550e-03   1.16e-02  1.489e-03     0.0      74      55
60  9.1103351e-03  1.3003737e-03   9.11e-03  1.124e-03     0.0      62      54
70  7.0428828e-03  8.7661595e-04   7.04e-03  8.631e-04     0.0      61      52
80  5.2766773e-03  1.2850278e-03   5.28e-03  6.801e-04     0.0      56      45
90  4.1388815e-03  8.8715473e-04   4.14e-03  5.256e-04     0.0      52      43
100  2.6822888e-03  3.2119283e-03   2.68e-03  4.857e-04    -0.3      44      34
110  2.0866055e-03  2.5426187e-03   2.09e-03  3.686e-04    -1.2      37      29
120  1.5623253e-03  3.0195558e-04   1.56e-03  2.006e-04     0.0      30      23
130  1.2843009e-03  1.3879354e-04   1.28e-03  1.609e-04     0.0      18      17
140  9.7991180e-04  1.3605821e-04   9.80e-04  1.299e-04     0.0      14      14
150  9.5062116e-04  7.3619749e-05   9.51e-04  1.284e-04     0.0      14      14
152  9.5037298e-04  5.7557443e-06   9.50e-04  1.242e-04     0.0      14      14  1.7049998e+01
160  1.1103540e-04  2.3479922e-05   1.11e-04  1.431e-05     0.0      14      14
164  9.5001274e-05  2.1214042e-05   9.50e-05  1.217e-05     0.0      14      14

EXIT -- Found a root

Products with A     :     266        Total time   (secs) :     0.2
Products with A^H   :     165        Project time (secs) :     0.1
Newton iterations   :       4        Mat-vec time (secs) :     0.0
Line search its     :     101        Subspace iterations :       0



## Pareto Frontier

Sample the Pareto frontier at 100 points:

$$phi(tau) = min. ||Ax-b||_2 \quad subject \quad to \quad ||x|| <= \tau$$


In [9]:

b = A.dot(x0)
x = np.zeros(n)
tau = np.linspace(0, 1.05 * np.linalg.norm(x0, 1), 100)
tau[0] = 1e-10
phi = np.zeros(tau.size)

for i in range(tau.size):
x,r,grad,info = spgl1(A, b, tau[i], 0, x, iter_lim=1000)
phi[i] = np.linalg.norm(r)

plt.figure()
plt.plot(tau,phi, '.')
plt.title('Pareto frontier')
plt.xlabel('||x||_1')
plt.ylabel('||Ax-b||_2');






## Weighted BP

Solve

$$min. ||y||_1 \quad subject \quad to \quad AW^{-1}y = b$$

and the weighted basis pursuit (BP) problem:

$$min. ||Wx||_1 \quad subject \quad to \quad Ax = b$$

followed by setting $y = Wx$.



In [10]:

# Sparsify vector x0 a bit more to get exact recovery
k = 9
x0 = np.zeros(n)
x0[p[0:k]] = np.random.randn(k)

# Set up weights w and vector b
w = np.random.rand(n) + 0.1 # Weights
b = A.dot(x0/w) # Signal

# Solution
x,resid,grad,info = spg_bp(A, b, **dict(iter_lim=1000, weights=w))

# Reconstructed solution, with weighting
x1 = x * w

plt.figure()
plt.plot(x1,'b')
plt.plot(x0,'ro')
plt.legend(('Coefficients','Original coefficients'))
plt.title('Weighted Basis Pursuit');






## MMV

Solve the multiple measurement vector (MMV) problem

$$(1) \quad min. ||Y||_{1,2} \quad subject \quad to \quad AW^{-1}Y = B$$

and the weighted MMV problem (weights on the rows of X):

$$(2) \quad min. ||WX||_{1,2} \quad subject \quad to \quad AX = B$$

followed by setting $Y = WX$.



In [11]:

# Create problem
m = 100
n = 150
k = 12
l = 6;
A = np.random.randn(m, n)
p = np.random.permutation(n)[:k]
X0 = np.zeros((n, l))
X0[p, :] = np.random.randn(k, l)

weights = 3 * np.random.rand(n) + 0.1
W = 1/weights * np.eye(n)

B = A.dot(W).dot(X0)

# Solve unweighted version
x_uw, _, _, _ = spg_mmv(A.dot(W), B, 0, **dict(verbosity=1))

# Solve weighted version
x_w, _, _, _ = spg_mmv(A, B, 0, **dict(verbosity=2, weights=weights))
x_w = spdiags(weights, 0, n, n).dot(x_w)

# Plot results
plt.figure()
plt.plot(x_uw[:, 0], 'b-', label='Coefficients (1)')
plt.plot(x_w[:, 0], 'g--', label='Coefficients (2)')
plt.plot(X0[:, 0], 'ro', label='Original coefficients')
plt.legend()
plt.title('Weighted Basis Pursuit with Multiple Measurement Vectors');

plt.figure()
plt.plot(x_uw[:, 1], 'b', label='Coefficients (1)')
plt.plot(x_w[:, 1], 'g--', label='Coefficients (2)')
plt.plot(X0[:, 1], 'ro', label='Original coefficients')
plt.legend()
plt.title('Weighted Basis Pursuit with Multiple Measurement Vectors');




================================================================================
SPGL1
================================================================================
No. rows              :      600
No. columns           :      900

Initial tau           : 0.00e+00
Two-norm of b         : 1.39e+02

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :     6000

Linesearch failed with error 1. Damping max BB scaling to 10000.0
Linesearch failed with error 1. Damping max BB scaling to 1000.0
Linesearch failed with error 1. Damping max BB scaling to 100.0
Linesearch failed with error 1. Damping max BB scaling to 10.0
Linesearch failed with error 1. Damping max BB scaling to 1.0
Linesearch failed with error 1. Damping max BB scaling to 0.1

EXIT -- Found a BP solution

Products with A     :    7205        Total time   (secs) :     9.4
Products with A^H   :    3762        Project time (secs) :     5.9
Newton iterations   :      44        Mat-vec time (secs) :     0.3
Line search its     :    3431        Subspace iterations :       0

================================================================================
SPGL1
================================================================================
No. rows              :      600
No. columns           :      900

Initial tau           : 0.00e+00
Two-norm of b         : 1.39e+02

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :     6000

iterr      Objective   Relative Gap  Rel Error      gnorm   stepg   nnz_x   nnz_g     tau

0  1.3889379e+02  0.0000000e+00   1.00e+00  5.510e+03     0.0       0       0  3.5012369e+00
1  1.0633151e+02  1.8286125e+00   1.00e+00  3.904e+03    -0.6     254       0
2  6.1884091e+01  6.9644415e-01   1.00e+00  1.352e+03     0.0      90       0
3  5.0592369e+01  1.9932604e-01   1.00e+00  4.709e+02     0.0      35       0
4  5.0412257e+01  1.3323811e-02   1.00e+00  3.475e+02     0.0      24       0
5  5.0411302e+01  1.8745545e-03   1.00e+00  3.463e+02     0.0      24       0  1.0839677e+01
6  2.9570102e+01  3.2810149e+00   1.00e+00  1.547e+02     0.0     228       0
7  2.7430718e+01  6.6709468e-01   1.00e+00  8.946e+01     0.0     208       2
8  2.7037539e+01  5.1896596e-01   1.00e+00  7.331e+01     0.0     151       1
9  2.6950560e+01  1.2789858e-01   1.00e+00  6.296e+01     0.0     147       4
10  2.6937796e+01  1.7817903e-01   1.00e+00  6.117e+01     0.0     143       4
12  2.6933994e+01  2.8948151e-02   1.00e+00  5.734e+01     0.0     142       4  2.3490402e+01
20  7.5332781e+00  1.3947307e+00   1.00e+00  1.075e+01     0.0     315      15
27  7.3695964e+00  1.1665351e-01   1.00e+00  9.453e+00     0.0     255      19  2.9235496e+01
30  2.1176118e+00  2.4175540e+02   1.00e+00  1.940e+01     0.0     679       0
40  7.6551591e-01  2.9799392e+02   7.66e-01  1.059e+01    -0.6     399       0
50  3.0701143e-01  1.1822538e+01   3.07e-01  6.696e-01     0.0     282       2
60  1.6413582e-01  1.9704894e+00   1.64e-01  2.711e-01     0.0      81       9
63  1.6390110e-01  1.5534069e+00   1.64e-01  2.572e-01     0.0      85       9  2.9339961e+01
70  4.4473829e-02  8.5733063e+00   4.45e-02  3.345e-01    -0.6      75       0
80  3.4576737e-02  1.4939973e-01   3.46e-02  4.737e-02     0.0      72      16
85  3.4543450e-02  5.9818837e-03   3.45e-02  4.277e-02     0.0      72      18  2.9367861e+01
90  7.9089320e-03  9.0402304e-01   7.91e-03  3.603e-02     0.0      72       0
100  2.4799153e-03  3.8387108e-01   2.48e-03  1.488e-02     0.0      72       0
110  1.3061292e-03  5.1282683e-01   1.31e-03  1.807e-02    -0.9      72       0
120  5.2941042e-04  1.7101981e-02   5.29e-04  9.866e-04     0.0      72      66
130  2.8342231e-04  9.9015771e-03   2.83e-04  5.550e-04     0.0      72      72
140  1.3217629e-04  1.9322320e-02   1.32e-04  7.644e-04     0.0      72      72
148  7.5109415e-05  5.3293672e-03   7.51e-05  2.415e-04     0.0      72      72

EXIT -- Found a BP solution

Products with A     :     192        Total time   (secs) :     0.2
Products with A^H   :     149        Project time (secs) :     0.1
Newton iterations   :       6        Mat-vec time (secs) :     0.0
Line search its     :      43        Subspace iterations :       0



## MMV with non-negative solution



In [12]:

# Create problem
m = 100
n = 150
k = 12
l = 6;
A = np.random.randn(m, n)
p = np.random.permutation(n)[:k]
X0 = np.zeros((n, l))
X0[p, :] = np.abs(np.random.randn(k, l))

B = A.dot(X0)

X, _, _, _ = spg_mmv(A, B, 0, iter_lim=10, verbosity=1)
XNN, _, _, _ = spg_mmv(A, B, 0, iter_lim=10, verbosity=1,
project=norm_l12nn_project,
primal_norm=norm_l12nn_primal,
dual_norm=norm_l12nn_dual)
print('Negative X:', np.any(X))
print('Negative XNN:', np.any(XNN))

# Plot results
plt.figure()
plt.plot(X[:, 0], 'b-', label='Coefficients')
plt.plot(XNN[:, 0], 'g--', label='Coefficients NN')
plt.plot(X0[:, 0], 'ro', label='Original coefficients')
plt.legend()
plt.title('Weighted Basis Pursuit with Multiple Measurement Vectors');

plt.figure()
plt.plot(X[:, 1], 'b', label='Coefficients')
plt.plot(XNN[:, 1], 'g--', label='Coefficients NN')
plt.plot(X0[:, 1], 'ro', label='Original coefficients')
plt.legend()
plt.title('Weighted Basis Pursuit with Multiple Measurement Vectors');




================================================================================
SPGL1
================================================================================
No. rows              :      600
No. columns           :      900

Initial tau           : 0.00e+00
Two-norm of b         : 8.07e+01

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :       10

ERROR EXIT -- Too many iterations

Products with A     :      13        Total time   (secs) :     0.0
Products with A^H   :      11        Project time (secs) :     0.0
Newton iterations   :       2        Mat-vec time (secs) :     0.0
Line search its     :       2        Subspace iterations :       0

================================================================================
SPGL1
================================================================================
No. rows              :      600
No. columns           :      900

Initial tau           : 0.00e+00
Two-norm of b         : 8.07e+01

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :       10

ERROR EXIT -- Too many iterations

Products with A     :      13        Total time   (secs) :     0.0
Products with A^H   :      11        Project time (secs) :     0.0
Newton iterations   :       2        Mat-vec time (secs) :     0.0
Line search its     :       2        Subspace iterations :       0
Negative X: True
Negative XNN: True



## LSQR

Let's finally try to compare the internal lsqr with scipy lsqr and perform sgpl1 with subspace minimization



In [13]:

def Aprodfun(A, x, mode):
if mode == 1:
y = np.dot(A,x)
else:
return np.dot(np.conj(A.T), x)
return y




In [14]:

n = 10
m = 20
A = np.random.normal(0, 1, (m, n))
Aprod = lambda x, mode: Aprodfun(A, x, mode)
x = np.ones(n)
y = A.dot(x)

damp       = 1e-5
aTol       = 1e-5
bTol       = 1e-5
conLim     = 1e12
itnMaxLSQR = 100
showLSQR   = 2

xinv, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var = \
lsqr(m, n, Aprod, y, damp, aTol, bTol, conLim, itnMaxLSQR, showLSQR)

xinv_sp, istop_sp, itn_sp, r1norm_sp, r2norm_sp, anorm_sp, acond_sp, arnorm_sp, xnorm_sp, var = \
splsqr(A, y, damp, aTol, bTol, conLim, itnMaxLSQR, showLSQR)

print('istop=%d, itn=%d, r1norm=%.2f, '
'r2norm=%.2f, anorm=%.2f, acond=%.2f, arnorm=%.2f, xnorm=%.2f' \
%(istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm))

print('istop=%d, itn=%d, r1norm=%.2f, '
'r2norm=%.2f, anorm=%.2f, acond=%.2f, arnorm=%.2f, xnorm=%.2f' \
%(istop_sp, itn_sp, r1norm_sp, r2norm_sp, anorm_sp, acond_sp, arnorm_sp, xnorm_sp))

plt.plot(x, lw=8)
plt.plot(xinv, '--g', lw=4)
plt.plot(xinv_sp, '--r')
plt.ylim(0, 2);




LSQR            Least-squares solution of  Ax = b
The matrix A has       20 rows  and       10 cols
damp = 1.00000000000000e-05    wantvar =        1
atol = 1.00e-05                 conlim = 1.00e+12
btol = 1.00e-05                 itnlim =      100

Itn      x(1)       r1norm     r2norm  Compatible   LS      Norm A   Cond A
0  0.00000e+00  8.439e+00  8.439e+00   1.0e+00  4.5e-01
1  4.29387e-01  4.857e+00  4.857e+00   5.8e-01  3.8e-01  4.6e+00  1.0e+00
2  9.26782e-01  3.020e+00  3.020e+00   3.6e-01  3.8e-01  5.8e+00  2.4e+00
3  1.06774e+00  1.992e+00  1.992e+00   2.4e-01  1.9e-01  8.6e+00  4.2e+00
4  1.09945e+00  9.175e-01  9.175e-01   1.1e-01  9.4e-02  9.6e+00  5.8e+00
5  9.92975e-01  4.452e-01  4.452e-01   5.3e-02  4.4e-02  1.0e+01  7.1e+00
6  1.01319e+00  1.885e-01  1.885e-01   2.2e-02  1.3e-02  1.1e+01  8.4e+00
7  1.01872e+00  1.020e-01  1.020e-01   1.2e-02  1.1e-02  1.1e+01  9.8e+00
8  1.01012e+00  4.644e-02  4.644e-02   5.5e-03  3.6e-03  1.2e+01  1.1e+01
9  9.99975e-01  1.085e-03  1.085e-03   1.3e-04  7.2e-05  1.3e+01  1.3e+01
10  1.00000e+00  1.499e-10  3.162e-05   3.7e-06  9.0e-13  1.3e+01  1.4e+01

LSQR finished
Ax - b is small enough, given atol, btol

istop =       1   r1norm = 1.5e-10   anorm = 1.3e+01   arnorm = 2.8e-11
itn   =      10   r2norm = 3.2e-05   acond = 1.4e+01   xnorm  = 3.2e+00

nprodA 10
nprodA 11

LSQR            Least-squares solution of  Ax = b
The matrix A has       20 rows  and       10 cols
damp = 1.00000000000000e-05   calc_var =        0
atol = 1.00e-05                 conlim = 1.00e+12
btol = 1.00e-05               iter_lim =      100

Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
0  0.00000e+00   8.439e+00  8.439e+00    1.0e+00  4.5e-01
1  4.29387e-01   4.857e+00  4.857e+00    5.8e-01  5.4e-01   4.6e+00  1.0e+00
2  9.26782e-01   3.020e+00  3.020e+00    3.6e-01  6.9e-01   5.8e+00  2.4e+00
3  1.06774e+00   1.992e+00  1.992e+00    2.4e-01  3.5e-01   8.6e+00  4.2e+00
4  1.09945e+00   9.175e-01  9.175e-01    1.1e-01  3.4e-01   9.6e+00  5.8e+00
5  9.92975e-01   4.452e-01  4.452e-01    5.3e-02  3.0e-01   1.0e+01  7.1e+00
6  1.01319e+00   1.885e-01  1.885e-01    2.2e-02  2.0e-01   1.1e+01  8.4e+00
7  1.01872e+00   1.020e-01  1.020e-01    1.2e-02  2.9e-01   1.1e+01  9.8e+00
8  1.01012e+00   4.644e-02  4.644e-02    5.5e-03  2.0e-01   1.2e+01  1.1e+01
9  9.99975e-01   1.085e-03  1.085e-03    1.3e-04  1.7e-01   1.3e+01  1.3e+01
10  1.00000e+00   1.498e-10  3.162e-05    3.7e-06  4.4e-08   1.3e+01  1.4e+01

LSQR finished
Ax - b is small enough, given atol, btol

istop =       1   r1norm = 1.5e-10   anorm = 1.3e+01   arnorm = 1.8e-11
itn   =      10   r2norm = 3.2e-05   acond = 1.4e+01   xnorm  = 3.2e+00

istop=1, itn=10, r1norm=0.00, r2norm=0.00, anorm=12.74, acond=14.39, arnorm=0.00, xnorm=3.16
istop=1, itn=10, r1norm=0.00, r2norm=0.00, anorm=12.74, acond=14.39, arnorm=0.00, xnorm=3.16



## Subspace minimization in SPGL1

And use subspace minimization in SPGL1



In [15]:

# Create random m-by-n encoding matrix and sparse vector
np.random.seed(0)

m = 50
n = 128
k = 14
[A, Rtmp] = np.linalg.qr(np.random.randn(n,m),'reduced')
A = A.T
p = np.random.permutation(n)
p = p[0:k]
x0 = np.zeros(n)
x0[p] = np.random.randn(k)

# Basis pursuit with subspace minimization
b = A.dot(x0) # signal
x,resid,grad,info = spg_bp(A, b, subspace_min=False, verbosity=2)
x,resid,grad,info_sub = spg_bp(A, b, subspace_min=True, verbosity=2)

plt.figure()
plt.plot(np.arange(info['niters']), info['rnorm2']/max(info['rnorm2']), '.-k',
label='without subspace min')
plt.plot(np.arange(info_sub['niters']), info_sub['rnorm2']/max(info_sub['rnorm2']), '.-r',
label='with subspace min')
plt.xlabel(r'#iter')
plt.ylabel(r'$||r||_2$')
plt.legend();




================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 2.20e+00

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :      500

iterr      Objective   Relative Gap  Rel Error      gnorm   stepg   nnz_x   nnz_g     tau

0  2.1958779e+00  0.0000000e+00   1.00e+00  8.769e-01     0.0       0       0  5.4988967e+00
1  1.3446156e+00  1.7302554e+00   1.00e+00  5.187e-01    -0.6      73       1
2  1.1528827e+00  9.9229525e-01   1.00e+00  3.836e-01     0.0      42       1
3  9.9709309e-01  3.2255337e-01   9.97e-01  2.181e-01     0.0      20       1
4  9.6583740e-01  2.7327432e-01   9.66e-01  2.254e-01     0.0      16       1
5  9.5588835e-01  2.1110135e-01   9.56e-01  1.911e-01     0.0      16       1
6  9.5213530e-01  2.5629885e-02   9.52e-01  1.659e-01     0.0      14       2
7  9.5182348e-01  1.7721386e-02   9.52e-01  1.637e-01     0.0      14       1
8  9.5159616e-01  9.3406208e-03   9.52e-01  1.622e-01     0.0      15       3
9  9.5145263e-01  1.0863959e-02   9.51e-01  1.618e-01     0.0      15       1
10  9.5156884e-01  3.0524527e-02   9.52e-01  1.663e-01     0.0      15       1
12  9.5144392e-01  1.9489737e-04   9.51e-01  1.601e-01     0.0      15      15  1.1153966e+01
20  1.6697472e-01  3.9246694e-02   1.67e-01  2.409e-02     0.0      55       2
30  1.4302223e-01  2.6788922e-02   1.43e-01  2.031e-02     0.0      49       2
40  1.1838948e-01  1.2786004e-02   1.18e-01  1.639e-02     0.0      46       6
50  1.1249910e-01  1.0412511e-02   1.12e-01  1.566e-02     0.0      37       7
60  1.1019625e-01  2.7113875e-03   1.10e-01  1.491e-02     0.0      35      35
62  1.1010028e-01  1.2995106e-02   1.10e-01  1.589e-02    -0.6      35       8  1.1917006e+01
70  1.6849338e-02  3.1194081e-03   1.68e-02  2.240e-03     0.0      62      62
80  1.5317457e-02  2.2139356e-03   1.53e-02  2.009e-03     0.0      55      55
90  1.4746957e-02  1.0477555e-03   1.47e-02  1.874e-03     0.0      52      52
99  1.4453429e-02  6.3972686e-03   1.45e-02  2.286e-03    -1.5      50      48  1.2008372e+01
100  4.3212271e-03  2.6249409e-03   4.32e-03  6.832e-04     0.0      57      57
110  3.5781423e-03  9.4808433e-04   3.58e-03  5.005e-04     0.0      56      56
120  3.4592825e-03  1.3929204e-03   3.46e-03  5.255e-04    -0.9      55      55
130  3.3442716e-03  1.6906443e-03   3.34e-03  5.350e-04    -0.6      55      55
140  3.1951155e-03  9.7620717e-04   3.20e-03  4.581e-04     0.0      52      52
150  3.1168251e-03  6.8760165e-04   3.12e-03  4.174e-04     0.0      48      48
160  3.0878221e-03  1.0734103e-04   3.09e-03  3.768e-04     0.0      48      48
170  3.0746380e-03  7.1288110e-05   3.07e-03  3.728e-04     0.0      48      48
177  3.0676839e-03  2.8380331e-04   3.07e-03  3.861e-04    -0.6      48      48  1.2032743e+01
180  3.6513571e-04  2.2617330e-04   3.65e-04  5.520e-05     0.0      49      49
190  2.5715605e-04  7.5652034e-05   2.57e-04  3.430e-05     0.0      49      49
200  2.3615363e-04  9.7525090e-05   2.36e-04  3.434e-05    -0.9      49      49
210  2.2754188e-04  8.6893137e-05   2.28e-04  3.254e-05    -1.2      49      49
215  2.2521623e-04  8.3088258e-05   2.25e-04  3.171e-05     0.0      49      49  1.2034343e+01
217  5.3849511e-05  3.8461025e-05   5.38e-05  8.578e-06     0.0      49      49

EXIT -- Found a root

Products with A     :     336        Total time   (secs) :     0.1
Products with A^H   :     218        Project time (secs) :     0.1
Newton iterations   :       6        Mat-vec time (secs) :     0.0
Line search its     :     118        Subspace iterations :       0

================================================================================
SPGL1
================================================================================
No. rows              :       50
No. columns           :      128

Initial tau           : 0.00e+00
Two-norm of b         : 2.20e+00

Optimality tol        : 1.00e-04
Target objective      : 0.00e+00

Basis pursuit tol     : 1.00e-06
Maximum iterations    :      500

iterr      Objective   Relative Gap  Rel Error      gnorm   stepg   nnz_x   nnz_g     tau

0  2.1958779e+00  0.0000000e+00   1.00e+00  8.769e-01     0.0       0       0  5.4988967e+00
1  1.3446156e+00  1.7302554e+00   1.00e+00  5.187e-01    -0.6      73       1
2  1.1489843e+00  9.8325846e-01   1.00e+00  3.821e-01     0.0      41       1                S 10
3  9.9519410e-01  3.3525857e-01   9.95e-01  2.186e-01     0.0      19       1
4  9.6624285e-01  2.7618469e-01   9.66e-01  2.266e-01     0.0      16       1
5  9.5515543e-01  1.7793784e-01   9.55e-01  1.866e-01     0.0      16       1
6  9.5214477e-01  2.6408557e-02   9.52e-01  1.660e-01     0.0      14       2
7  9.5183833e-01  1.8398733e-02   9.52e-01  1.639e-01     0.0      14       1
8  9.5158096e-01  8.8651909e-03   9.52e-01  1.621e-01     0.0      15       3
9  9.5145458e-01  1.0671337e-02   9.51e-01  1.617e-01     0.0      15       3
10  9.5158443e-01  3.0607887e-02   9.52e-01  1.665e-01     0.0      15       1
12  9.5144389e-01  8.6736938e-05   9.51e-01  1.600e-01     0.0      15      15  1.1155030e+01 S  2
20  1.6372851e-01  3.9315399e-02   1.64e-01  2.365e-02     0.0      54       2
30  1.2650967e-01  2.2251007e-02   1.27e-01  1.818e-02     0.0      46       3
40  1.1177426e-01  1.1040318e-02   1.12e-01  1.584e-02     0.0      38       9
50  1.1038074e-01  3.2093537e-02   1.10e-01  1.781e-02    -0.9      35       3
60  1.0980987e-01  4.9109500e-04   1.10e-01  1.468e-02     0.0      34      34
61  1.0980865e-01  9.3021046e-04   1.10e-01  1.472e-02     0.0      34      34  1.1974405e+01 S  2
70  1.0297557e-02  2.2054106e-03   1.03e-02  1.348e-03     0.0      66      66
80  8.5678149e-03  2.9152140e-03   8.57e-03  1.237e-03     0.0      56      56
90  8.0017968e-03  4.6253015e-03   8.00e-03  1.313e-03     0.0      54      54
100  7.5630334e-03  3.7270005e-03   7.56e-03  1.205e-03    -0.9      52      52
110  7.2362447e-03  7.5933779e-04   7.24e-03  9.419e-04     0.0      50      50
115  7.2201645e-03  5.0373057e-04   7.22e-03  9.168e-04    -0.6      50      50  1.2031267e+01
120  7.8540090e-04  1.1638891e-03   7.85e-04  1.676e-04    -0.6      53      53
130  5.2568119e-04  7.4888173e-05   5.26e-04  6.613e-05     0.0      52      52
140  4.7045197e-04  9.7612646e-05   4.70e-04  6.250e-05    -0.9      51      51
150  4.3029018e-04  6.4297913e-05   4.30e-04  5.644e-05     0.0      52      52
160  4.2071938e-04  8.7009640e-05   4.21e-04  5.640e-05    -1.2      52      52
170  4.1432024e-04  8.9704424e-05   4.14e-04  5.566e-05     0.0      52      52
180  4.0937269e-04  1.8428929e-05   4.09e-04  4.942e-05     0.0      52      52
190  4.0579936e-04  1.1508572e-04   4.06e-04  5.735e-05    -1.2      52      52
200  4.0230922e-04  2.9340751e-05   4.02e-04  4.924e-05     0.0      52      52
210  3.9695255e-04  6.4319323e-05   3.97e-04  5.114e-05    -0.6      52      52
220  3.9414669e-04  1.3021692e-05   3.94e-04  4.693e-05     0.0      52      52
222  3.9391597e-04  3.8706367e-05   3.94e-04  4.895e-05    -0.9      52      52  1.2034437e+01
223  5.1886239e-05  7.9437681e-05   5.19e-05  1.074e-05     0.0      52      52                S  8

EXIT -- Found a root

Products with A     :     583        Total time   (secs) :     0.2
Products with A^H   :     647        Project time (secs) :     0.1
Newton iterations   :       5        Mat-vec time (secs) :     0.0
Line search its     :     159        Subspace iterations :       8