# Compressed Sensing using StructuredOptimization.jl

**URL:** <https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098>\
**Category:** Optimization (Mathematical)\
**Created:** [February 25, 2020, 12:48am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098 "2020-02-25T00:48:14Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [February 25, 2020, 12:48am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/1 "2020-02-25T00:48:14Z")

</div>

Here is something cool I did: [gist](https://gist.github.com/lstagner/79c1feb91210b8829c00552ce0b392a8)

```julia
using LinearAlgebra
using Random
using FFTW
using PyPlot
using StructuredOptimization

N = 256
P = 5
K = 32

# Pick out P random frequencies
freq = randperm(Int(N/2)).-1
freq = freq[1:P]

# Construct Signal
n = 0:N-1
x = zeros(N)
for f in freq
    x .+= sin.((2pi*f/N)*n)
end

# Construct DFT Matrix
Psi = hcat(fft.(eachcol(I(N)))...)
Psi_inv = conj(Psi)/N
X = Psi*x # Fourier Components

# Pick K Random Data
x_m = zeros(N);
q = randperm(N);
q = q[1:K];
x_m[q] .= x[q]

#Construct system of linear equations
A = Psi_inv[q, :]; 
y = A*X; # y === x_m[q]

# Use StructuredOptimization.jl to optimize the system with a L1 norm
xx = Variable(Complex{Float64}, N)
lambda = 1e-3*norm(A'*y,Inf)
@minimize ls(A*xx - y) + lambda*norm(xx,1)

#Plot it
fig, ax = plt.subplots(ncols=2)
fig.set_size_inches(12,4)
ax[1].plot(1:N,x,label="True Signal")
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Amplitude")
ax[1].plot((1:N)[q], x_m[q],"x",label="Data Points")
ax[1].plot(1:N, Psi_inv*(~xx),ls="--",label="Reconstructed Signal")
ax[1].legend()

ax[2].plot(abs.(X),label ="True Frequency")
ax[2].plot(abs.(~xx),ls="--",color="green", label="Reconstructed Frequency")
ax[2].set_xlim(1,N/2)
ax[2].set_xlabel("Frequency")
ax[2].set_ylabel("Amplitude")
ax[2].legend(loc=2)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/7/470754a9bf1863b4cd9815ab764b471e6ab34642.png)

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [February 25, 2020, 1:31am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/2 "2020-02-25T01:31:09Z")

</div>

I would love to have a blog of compressed sensing examples! I could probably help contribute 1 or 2. Have you seen [GitHub - JuliaFirstOrder/ProximalOperators.jl: Proximal operators for nonsmooth optimization in Julia](https://github.com/kul-forbes/ProximalOperators.jl) ?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 25, 2020, 4:49am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/3 "2020-02-25T04:49:49Z")

</div>

I have implemented some sparse spectral estimation methods here

> **[GitHub - baggepinnen/LPVSpectral.jl: Least-squares (sparse) spectral...](https://github.com/baggepinnen/LPVSpectral.jl)**
>
> Least-squares (sparse) spectral estimation and (sparse) LPV spectral decomposition. - GitHub - baggepinnen/LPVSpectral.jl: Least-squares (sparse) spectral estimation and (sparse) LPV spectral decom...

---

<div class="post-metadata">

**Author:** ![tobias.knopp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobias.knopp/32/7551_2.png) [@tobias.knopp](https://discourse.julialang.org/u/tobias.knopp)\
**Post date:** [February 25, 2020, 11:11am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/4 "2020-02-25T11:11:11Z")

</div>

We are working on compressed sensing in the field of image reconstruction. This package contains the basic solver [GitHub - tknopp/RegularizedLeastSquares.jl](https://github.com/tknopp/RegularizedLeastSquares.jl) where one can plug in different imaging operators and sparsity constraints. This is used in our MRI package [GitHub - MagneticResonanceImaging/MRIReco.jl: Julia Package for MRI Reconstruction](https://github.com/MagneticResonanceImaging/MRIReco.jl)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 25, 2020, 11:22am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/5 "2020-02-25T11:22:14Z")

</div>

That looks like an awesome package!

---

<div class="post-metadata">

**Author:** ![nantonel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nantonel/32/2889_2.png) [@nantonel](https://discourse.julialang.org/u/nantonel)\
**Post date:** [February 25, 2020, 11:25am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/6 "2020-02-25T11:25:45Z")

</div>

There are a few [demos](https://kul-forbes.github.io/StructuredOptimization.jl/stable/demos/) of [StructuredOptimization](https://kul-forbes.github.io/StructuredOptimization.jl) doing similar things! Most of them deal with compressed sensing.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 25, 2020, 12:50pm UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/7 "2020-02-25T12:50:39Z")

</div>

I’ve actually played around with StructuredOptimization today, I’m solving a robust PCA problem, for which I have a custom solver that works quite well. I would like to experiment with variations of this problem, but am not getting StructuredOptimization to converge unless I have a really small problem size. The sizes below is borderline working. Maybe you have some idea on how to choose the solver and params for this kind of problem?

```julia
using TotalLeastSquares, Random, DSP

## Generate chirp signal
T = 2^14
e = randn(T)
fs = 8_000
t = range(0,step=1/fs, length=T)
chirp_f = LinRange(2000, 2200, T)
y = sin.(2pi .* chirp_f .* t )
yn = y + e
H = hankel(yn, 100)
@time TotalLeastSquares.rpca(H, tol=1e-3)
# about 6 seconds

using StructuredOptimization
λ = 1/sqrt(maximum(H))
L = Variable(size(H)...)
S = Variable(size(H)...)
@time @minimize ls(L+S-H) + λ*norm(S,1) st rank(L)<=4 with PANOC(tol=1e-2)
# > 1 minute with 60GiB allocations

```

---

<div class="post-metadata">

**Author:** ![nantonel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nantonel/32/2889_2.png) [@nantonel](https://discourse.julialang.org/u/nantonel)\
**Post date:** [February 25, 2020, 1:54pm UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/8 "2020-02-25T13:54:54Z")

</div>

Concerning speed, `rank` takes a truncated SVD. This is not done in place at the [moment](https://github.com/kul-forbes/ProximalOperators.jl/blob/master/src/functions/indBallRank.jl). I’m not sure if this is the bottleneck compared to TotalLeastSquares which seems to be exploiting the matrix Hankel structure as well. But surely the SVD is expensive since is taken at each iteration.

One thing you could try to avoid svd is to solve this instead:

```julia
using StructuredOptimization
n = 10
H = randn(n,n)
r = 4
U = Variable(randn(n,r))
Q = Variable(randn(r,n))
S = Variable(n,n)
lambda = 1/sqrt(maximum(H))
@minimize ls(U*Q+S-H) + lambda*norm(S,1)

```

Of course you can add some constraints/regularizations on `U` and `Q`.

About parameter tuning, I’m not that experienced about this problem. I’ve only worked on it with the demo and hand-tuned. For sure if you choose a rank that is too low and lambda too large you won’t be able to properly minimize the least squares term.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 25, 2020, 11:35pm UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/9 "2020-02-25T23:35:26Z")

</div>

`rpca` used in the example uses the `NuclearNorm` proximal operator so it also performs an SVD in each step. Is it possible to use `rank(L)` as a reglarization term or it can only be used as constraint in StructuredOptimization?

```julia
@time @minimize ls(L+S-H) + λ*norm(S,1) + rank(L)

```

throws an error about no prox available for `Rank`.

---

<div class="post-metadata">

**Author:** ![nantonel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nantonel/32/2889_2.png) [@nantonel](https://discourse.julialang.org/u/nantonel)\
**Post date:** [February 26, 2020, 11:39am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/10 "2020-02-26T11:39:10Z")

</div>

It is possible, the notation is the following:

```julia
using StructuredOptimization
n = 10
H = randn(n,n)
lambda1, lambdan = 1/sqrt(maximum(H)), 1e-3
L, S = Variable(size(H)...), Variable(size(H)...)
@minimize ls(L+S-H) + lambda1*norm(S,1) + lambdan*norm(L,*)

```

---

<div class="post-metadata">

**Author:** ![lostella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostella/32/356_2.png) [@lostella](https://discourse.julialang.org/u/lostella)\
**Post date:** [February 26, 2020, 12:04pm UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/11 "2020-02-26T12:04:44Z")

</div>

@baggepinnen for clarity, note that `norm(L, *)` will result in a nuclear norm penalty, not rank (which is maybe what you originally wanted)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 26, 2020, 12:42pm UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/12 "2020-02-26T12:42:06Z")

</div>

Sorry I was not being entirely clear, I was thinking about the nuclear norm of course 😃

---

<div class="post-metadata">

**Author:** ![marknzed](https://avatars.discourse-cdn.com/v4/letter/m/c5a1d2/32.png) [@marknzed](https://discourse.julialang.org/u/marknzed)\
**Post date:** [September 26, 2021, 11:35am UTC](https://discourse.julialang.org/t/compressed-sensing-using-structuredoptimization-jl/35098/13 "2021-09-26T11:35:59Z")

</div>

Hi,

Sometimes (when re-evaluating q = randperm(N)) the reconstructed signal has a very good match and sometimes it is like the attached. Is this expected?

![Screen Shot 2021-09-26 at 13.31.34](https://global.discourse-cdn.com/julialang/original/3X/e/b/eb720025c5638b664d1306fa3cb78e43f3ada9c9.png)