Our goal is to find out the solution with a least possible non-zero entry of all of the solutions. That is, the solution should give us as few non-zeros as possible. Are you wondering where this can be applied? There are plenty of applications for it. The areas where it can be applied are as follows:
- Signal representation
- Medical imaging
- Sparse channel estimation
Let's say that we have got a time signal. This signal is highly sparse, but we know a little bit about it as it has a few frequencies. Can you sense what it is from the earlier equation? Yes, it can be deemed as X.
Let's call this unknown signal X. Now, even though we don't know the whole signal, we can still make observations about it, or samples, as shown in the following code:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.optimize as spopt
import scipy.fftpack as spfft
import scipy.ndimage as spimg
import cvxpy as cvx
This will form a random equation:
x = np.sort(np.random.uniform(0, 15, 30))
y = 5 + 0.5 * x + 0.1 * np.random.randn(len(x))
Now, we need to fit the l1 norm. We get the following output:
l1 = lambda x0, x, y: np.sum(np.abs(x0[0] * x + x0[1] - y))
opt1 = spopt.fmin(func=l1, x0=[1, 1], args=(x, y))
Then, we need to fit the l2 norm. We get the following output:
l2 = lambda x0, x, y: np.sum(np.power(x0[0] * x + x0[1] - y, 2))
opt2 = spopt.fmin(func=l2, x0=[1, 1], args=(x, y))
y2 = y.copy()
y2[3] += 5
y2[13] -= 10
xopt12 = spopt.fmin(func=l1, x0=[1, 1], args=(x, y2))
xopt22 = spopt.fmin(func=l2, x0=[1, 1], args=(x, y2))
By summing up the two sinusoids, we get the following output:
n = 10000
t = np.linspace(0, 1/5, n)
y = np.sin(1250 * np.pi * t) + np.sin(3000 * np.pi * t)
yt = spfft.dct(y, norm='ortho')
plt.figure(figsize=[10,5])
plt.plot(t,y)
plt.title('Original signal')
plt.xlabel('Time (s)')
plt.ylabel('y')
Now, let's take the sample out of n:
m = 1000 # 10% sample
ran = np.random.choice(n, m, replace=False) # random sample of indices
t2 = t[ran]
y2 = y[ran]
Let's create the idct matrix operator:
# create idct matrix operator
A = spfft.idct(np.identity(n), norm='ortho', axis=0)
A = A[ran]
# do L1 optimization
vx = cvx.Variable(n)
objective = cvx.Minimize(cvx.norm(vx, 1))
constraints = [A*vx == y2]
prob = cvx.Problem(objective, constraints)
result = prob.solve(verbose=True)
The output for this is as follows:
To reconstruct the signal, we must do the following:
x = np.array(vx.value)
x = np.squeeze(x)
signal = spfft.idct(x, norm='ortho', axis=0)
That is how we reconstruct the signal.