Test the use of Forward-backward-like splitting for the resolution of a compressed sensing regularization
In [1]:
from __future__ import division
import time
import numpy as np
import scipy.linalg as lin
import pylab as pl
from locore.algorithms import forward_backward
from locore.operators import soft_thresholding
Parameters
In [2]:
n = 600
p = n // 4
la = 1.0 # regularization parameter
Matrix and observations
In [3]:
A = np.random.randn(p, n)
y = np.random.randn(p, 1)
List of benchmarked algorithms
In [4]:
methods = ['fb', 'fista', 'nesterov']
Operator callbacks
In [5]:
F = lambda x: la * lin.norm(x, 1)
G = lambda x: 1 / 2 * lin.norm(y - np.dot(A, x)) ** 2
prox_f = lambda x, tau: soft_thresholding(x, la * tau)
grad_g = lambda x: np.dot(A.T, np.dot(A, x) - y)
Run
In [6]:
L = lin.norm(A, 2) ** 2 # Lipschitz constant
callback = lambda x: F(x) + G(x)
maxiter = 1000
res = np.zeros((maxiter, len(methods)))
i = 0
for method in methods:
t1 = time.time()
x, fx = forward_backward(prox_f, grad_g, np.zeros((n, 1)), L,
maxiter=maxiter, method=method,
full_output=1, retall=0, callback=callback)
t2 = time.time()
print "[" + method + "]: Performed 1000 iterations in " \
+ str(t2 - t1) + "seconds."
res[:, i] = fx
i += 1
Show
In [7]:
e = np.min(res.flatten())
pl.loglog(res[:(maxiter // 10), :] - e)
pl.legend(methods)
pl.grid(True, which="both", ls="-")
pl.tight_layout()
pl.show()