Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am playing with some Compressed Sensing (think single pixel camera) applications and would like to have a Mathematica equivalent of a Matlab package call Convex Optimization (CVX). Specifically, I am solving an underdetermined, constrained, L1Norm minimization problem. I have been using NMinimize, but it very slow compared to the Matlab code written by a colleague. I don’t want to give him the pleasure of thinking Matlab is superior. I am sure that the long run times are due to my poor programing skills. The solution is sparse but I do not know how to add that as an additional constraint. Anyway, seems I can’t find a good reference for compressed sensing in Mathematica. Any suggestions would be great!

As an example two data sets

1D

coef = Transpose[RandomReal[#, 7] & /@ {1, 10 Pi, 2 Pi}];
f[t_] = Total[(#[[1]] Sin[#[[2]] t + #[[3]]]) & /@ coef];
data1d = Table[f[t], {t, 0, 10, 0.01}];

Sample obtained using

RandomSample[data1d,n1]

2D

data2d = ImageData[ExampleData[{"TestImage", "Boat"}]];

1 pixel random sampling

Total[RandomSample[#, 256] & /@ data2d, 2]

I hope the 200 reputation points bounty don't go to waste.

share|improve this question
6  
For a convex optimization problem FindMinimum is likely to be much faster than NMinimize. –  Daniel Lichtblau Jul 31 '14 at 23:33
1  
@rhermans I know you won't get pinged, but since you seem interested in this, consider editing the question to add your own example signal and re-sample. While I probably won't be able to solve the problem, others might be more tempted to try if they had an example that you chose for its properties. I think such an edit would be within the SE guidelines for editing others' questions. –  Michael E2 Jan 28 at 13:59
1  
FWIW, the Matlab CVX package is open-source: cvxr.com/cvx –  Michael E2 Feb 2 at 13:50
1  
Can you post, what you already tried (the NMinimize-approach) and some timing comparisons, so as to know what we are up against? –  Jinxed Feb 2 at 17:51
2  
@rhermans I assume the examples are yours. If so, I cannot say anything unless you post some code, as I have no idea what it is you are attempting to optimize. (If you don't wish to edit someone else's question, could post it as a new one. Done correctly, it would have sufficiently more info that it would not be regarded as a duplicate.) –  Daniel Lichtblau Feb 2 at 18:07

1 Answer 1

I am playing with some Compressed Sensing (think single pixel camera) applications and would like to have a Mathematica equivalent of a Matlab package call Convex Optimization (CVX). [...] very slow compared to the Matlab code written by a colleague. I don’t want to give him the pleasure of thinking Matlab is superior

CVX is the result of years of theoretical and applied research, a book on convex optimization and a company focused on researching, developing and supporting convex optimization tools. You simply cannot create a Mathematica clone overnight that parallels the performance and features of CVX, and certainly not via a question on Stack Exchange! :)

There are plenty (way more than you'll ever need!) of examples with code for doing compressed sensing/L1 optimizations in MATLAB for different constraints and your best bet would be to leverage those existing scripts.

Use CVX via MATLink

The best way to do convex optimization (beyond what NMinimize and friends allow you) in Mathematica if you also have a MATLAB license is to use CVX via MATLink. In fact, I co-wrote it with Szabolcs specifically because I personally needed to use CVX + Mathematica.

As an example workflow of using CVX with Mathematica, we can solve this convex optimization problem as follows:

MEvaluate["
    n = 6;
    m = 40;
    randn('state',0);

    % generate 50 ponts ui, vi
    u = linspace(-1,1,m);
    v = 1./(5+40*u.^2) + 0.1*u.^3 + 0.01*randn(1,m);

    A = vander(u');
    A = A(:,m-n+[1:n]);     

    % L-infty fit
    cvx_begin quiet
       variable x_inf(n)
       minimize (norm(A*x_inf - v', inf))
    cvx_end
"]

CVX is quite verbose and it is a bit more efficient to use quiet when running it via MATLink. There are other CVX global variables (starting with cvx_) that you can use to programmatically check which solver was used, whether the problem was solved or if it is infeasible, etc. rather than relying on the textual output. If you are doing active algorithm development using CVX, then it might be better to use MATLAB directly for this part and then use MATLink once you are satisfied and have frozen the routine.

Next, we transfer the optimal solution and the original quantities to Mathematica using MGet:

{u, v, A, xinf} = MGet[{"u", "v", "A", "x_inf"}];

Finally, we compute the optimal $L_2$ norm solution in Mathematica, construct the polynomial using Horner's method and plot them in Mathematica:

With[{hornerPoly = Function[{var, x}, Fold[# var + #2 &, x]]},
    Show[
        ListPlot[Transpose@{u, v}, PlotStyle -> {AbsolutePointSize[5], Gray}],
        Plot[{hornerPoly[x, x2], hornerPoly[x, xinf]}, {x, -1, 1}, Evaluated -> True,
            PlotLegends -> {"!(*SubscriptBox[(L), (2)])",  "!(*SubscriptBox[(L), ([Infinity])])"}]
    ]
]

Compressed sensing/Sparse representations in Mathematica

If you don't have access to MATLAB, then you'll have to implement the algorithms yourself. You'll find NMinimize thoroughly lacking when you want to implement $\ell_1$ minimization and enforce sparsity or do $\ell_p$ minimizations for $p < 1$, etc. There's no magic way around actually learning the theory behind sparse representations yourself, learning the standard algorithms which are the building blocks to developing more sophisticated algorithms for your specific problem, understanding standard computational tricks to enforce sparsity, restricted isometry property, etc. I can recommend the following books which I found quite helpful:

Not to leave you completely in the dark, here are my implementations of a few basic pursuit algorithms (written for version 9). Also note that I did not write this (2-3 years old at this point) with the intent of sharing with others or as an expository piece of code, so it might be quite cryptic but reasonably idiomatic. Nevertheless, they're based on the pseudocode in Chapter 3 of Elad's book, so if you're using that, you might be able to follow.

Matching Pursuit

MatchingPursuit[A_?MatrixQ, b_, epsilon_] := Module[
    {solution, residual},

    solution[] = ConstantArray[0, Last@Dimensions@A];
    solution[z_] := With[
        {pos = First@Position[z, Max@z]},
        solution[] = ReplacePart[solution[], pos -> First@(solution[][[ pos ]] + 
            Transpose[A[[;;, pos]]].residual[])];
        pos
    ];

    residual[] = b;
    residual[pos_] := residual[] = residual[] - (A[[;;, pos]].Transpose@A[[;;, pos]]).residual[];

    Sow@NestWhile[Composition[residual,solution,Abs[Transpose@A.#]&], b, Flatten@#.Flatten@# > epsilon&];
    solution[] // Chop
]

Iteratively Re-weighted Least Squares (IRLS)

IRLS[A_?MatrixQ, b_, epsilon_, p_:1] := NestWhile[
    Composition[
        DiagonalMatrix@SparseArray@Flatten@Abs@#^p&, 
        #.Transpose@A.PseudoInverse[A.#.Transpose@A].b&, 
        #.#&],
    Last@Dimensions@A ~IdentityMatrix~ SparseArray,
    Sow@Norm[#-#2] > epsilon&, 2
]["NonzeroValues"] // Chop

Orthogonal Matching Pursuit

OrthogonalMatchingPursuit[A_?MatrixQ, b_, epsilon_] := Module[
    {support,residual},

    support[] = {};
    support[z_] := support[] = Flatten@Sort[support[] ~Join~ First@Position[z, Max@z]];

    residual := With[{AS = A[[;;,#]]}, b - (AS.PseudoInverse@AS).b]&;

    Sow@NestWhile[Composition[residual,support,Abs[Transpose@A.#]&], b, Flatten@#.Flatten@# > epsilon&];
    SparseArray[Thread[# -> Flatten[PseudoInverse@A[[;; , #]].b] &@support[]], Last@Dimensions@A] // Normal // Chop
]

Finally, here's a little snippet showing that each of these algorithms recover the sparse solution (again, these might make more sense if you're following Elad's book, but nevertheless):

SparseRandomVector[n_Integer,k_Integer] /; k < n := RandomSample[Table[RandomReal[], {k}] ~Join~ ConstantArray[0, n-k], n]
RandomMatrix[n_Integer,m_Integer] := RandomReal[NormalDistribution[], {n,m}]
NormalizeColumns[mat_?MatrixQ] := Normalize /@ Transpose[mat] // Transpose

SeedRandom[1234];
{#, MatrixForm@With[{A = NormalizeColumns@RandomMatrix[20, 40], 
       x = SparseRandomVector[40, 2], tol = 10^-4}, {x, #[A, A.x, tol]}]
 } & /@ {MatchingPursuit, OrthogonalMatchingPursuit, IRLS} // TableForm

share|improve this answer
    
@rhermans I missed your 200 bounty (was on vacation) but I hope the answer is helpful :) –  The Toad Feb 8 at 17:09
1  
I'm a bit confused by the remark on lacking l_1 optimization. That's linear programming. Do you mean handling of linear constraints with nonlinear objective (or vice versa)? –  Daniel Lichtblau Feb 8 at 20:26
    
@DanielLichtblau The remark included the clause following it — "and enforce sparsity". NMinimize handles $\ell_1$ minimization fine. In compressive sensing/sparse theory, one wishes to solve an $\ell_0$ minimization problem, but because that is NP-hard, you do $\ell_p,\ 0<p\le 1$ minimization under certain assumptions/conditions and with other possibly non-linear constraints. I guess my phrasing wasn't the best way to put it but nothing better comes to mind now... will update it if I think of something better. –  The Toad Feb 10 at 17:38
    
Okay. I was aware that l_1 was being used as a surrogate for l_0 for the problem at hand. It was the "with other possibly non-linear constraints" part I was missing. Makes sense now. –  Daniel Lichtblau Feb 10 at 17:41

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.