4
$\begingroup$

Hi All, I posted this on the math.se site, but this may be a better location.

I need a method of finding the maximum of a real valued trigonometric polynomial where I can trade accuracy for speed. The accepted answer to this question:

https://mathoverflow.net/questions/35538/the-maximum-of-a-real-trigonometric-polynomial

gives a method using semidefinite programming:

Let $f(x)=F(e^{ix})$ where $F(z)=\sum_{n=-N}^{N}c_n z^n$, with $c_n=\tfrac{1}{2}(a_n−i b_n)$ and $c_{−n}=\bar{c}_n$. Then $\min_x \,f(x)$ is equal to $c_0$ minus the value of the following semidefinite program: $\min_F tr(F)$ such that $F⪰0$, and $\sum_{p=k}^{N} F_{p,p−k}=c_k$ for $k=1,…,N$.

However, I don't understand what this means. Would anyone be able to explain in simpler terms and give an idea of how one would code this in MATLAB?

Update:

Some clarifications:

  • How does one implement a semidefinite program? Broad question I know! Perhaps better is what MATLAB package is best suited to this kind of semidefinite program problem? $N$ will be between 3 and 15 so I would like a quick solver for small matrices.
  • In the problem $F(z)$ is a function and $F$ is a matrix (I assume they are different)? When the SDP solves for $F$, how do I get the value of $z$ and $F(z)$ that corresponds to the maximum of $F(z)$?
Improve this question
$\endgroup$
7
  • $\begingroup$ geometrikal, welcome to SciComp! I agree with Brian here and encourage your to clarify your question. What aspects of the SDP would you like clarified? $\endgroup$
    – Geoff Oxberry
    Mar 19, 2013 at 16:29
  • $\begingroup$ @GeoffOxberry thanks for the welcome! My question shows my unfamiliarity with SDP i'm afraid. I've tried to clarify the question, but if some of the questions are too broad or ill-conceived some pointers to a good textbook or online resource would be welcome instead. $\endgroup$
    – geometrikal
    Mar 20, 2013 at 0:10
  • $\begingroup$ No problem! A couple other pointers: it's easier to answer a question that poses one (or two) questions. Your first bullet point, I think, is enough information for people to answer your question. Your remaining bullet points and the final sentence could be broken up into separate questions. The last two bullet points are probably one question, and that question could either fit here or at Math.SE, because that question doesn't require computation (your link is to MathOverflow, which is a different site). $\endgroup$
    – Geoff Oxberry
    Mar 20, 2013 at 4:57
  • $\begingroup$ By asking many questions at once, you risk answers that answer one well and others badly. (Disclaimer: this is purely subjective, and based on my experience on this site over the past year. It is your question, and I encourage you to do what you think is best.) $\endgroup$
    – Geoff Oxberry
    Mar 20, 2013 at 4:57
  • $\begingroup$ Hi @GeoffOxberry, I've updated the question. Thanks for the feedback. $\endgroup$
    – geometrikal
    Mar 20, 2013 at 5:03

3 Answers 3

Reset to default
4
$\begingroup$

I would not recommend using SDP to solve this problem. If roots is already obtaining the correct answer, I think you should stick with it or find some way to accelerate that approach; because the SDP method will be slower and less accurate. But just in case I am wrong about this, or just in case you insist on trying for other reasons, let me continue.

Johan's package, YALMIP/SDPT3 is certainly a good choice for this semidefinite program. So is CVX, which (full-disclosure) I co-authored. I don't think Johan got the constraints quite right, but the error is minor and I'm sure that will get fixed in due course. (EDIT: it has been fixed.) In the meanwhile I will share the CVX model.

Let's assume that you have stored the positive coefficients of $c$ in a length-$N+1$ MATLAB vector c, so that $c_k$ is in c(k+1) (because MATLAB indexes from 1, not 0). Then assuming the theorem you've cited is correct, the following CVX model should solve it:

cvx_begin sdp
    variable F(N+1,N+1) Hermitian
    maximize(c(1)-trace(F))
    for k = 1 : N,
        sum(diag(F,k)) == c(k+1)
    end
    F >= 0
cvx_end

Note that the objective function includes the $c_0$ term, and is subtracting the trace of $F$, so that step is included.

A disclaimer: while I'm sure I'm accurately representing the model you have stated in your question, I'm not 100% convinced that the model itself is correct. It's close, for sure---I'm familiar with the papers in question---but you would want to double-check if you go down this route.

EDITED TO ADD: Readers interested in potential applications for this type of model may wish to consult this recent paper, "Towards a Mathematical Theory of Super-Resolution", by Candès and Fernandez-Granda. They describe the use of trigonometric polynomial optimization to solve continuous-frequency sparse recovery applications.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for the answer. That paper sounds very interesting... Am I correct in my brief skim that it can be used for sparse recovery of the signal from its Fourier series coefficients? i.e. you don't need a vector of the discrete signal? $\endgroup$
    – geometrikal
    Mar 21, 2013 at 3:54
  • $\begingroup$ Ill check the answer against MATLAB roots. The reason I need speed is because I have to apply this at every pixel in an image. If there is a faster method that takes an initial guess that would be helpful too. Right now using roots takes about 2.3s for a 256x256 image using 4 cores $\endgroup$
    – geometrikal
    Mar 21, 2013 at 4:01
  • $\begingroup$ P.S. really great that there is a SE site like this. $\endgroup$
    – geometrikal
    Mar 21, 2013 at 4:02
  • $\begingroup$ What's a typical value of N? I am absolutely certain that CVX will not be fast enough---the modeling overhead alone will consume much more than 2.3s in total over the entire image. There's a cost to be paid to make the model look so pretty in MATLAB :-) Even if you bypassed CVX and called the SDP solver directly (e.g., SeDuMi, CSDP, SDPT3), I don't think there's one around that could do it that fast. You're talking 35 microseconds per pixel. $\endgroup$
    – Michael Grant
    Mar 21, 2013 at 12:10
  • $\begingroup$ N is 5 to 15 typically, although could be higher but then not every pixel will be evaluated. It's ok though - tried out the CVX code for the SDP in that paper and it solves another problem I've been having. :D Thanks again for the link. $\endgroup$
    – geometrikal
    Mar 21, 2013 at 13:12
3
$\begingroup$

If you install YALMIP (a modelling layer I develop) and a semidefinite solver (for instance SDPT3), the code would be, if I interpreted it correctly.

F = sdpvar(N+1,N+1,'hermitian','complex')
Constraints = F >= 0;
for k = 1:N
     Constraints = [Constraints, sum(diag(F,k)) == c(k+1)];
end
solvesdp(Constraints,trace(F));
double(F)

Note the shift in indices in F. It looks like the matrix F uses 0-based indexing which clashes with MATLABs 1-based.

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Main.HomePage

http://www.math.nus.edu.sg/~mattohkc/sdpt3.html

Improve this answer
$\endgroup$
7
  • $\begingroup$ Johan, you actually have to sum the $k$-th diagonal in each constraint. You've created N^2 constraints instead. $\endgroup$
    – Michael Grant
    Mar 20, 2013 at 15:45
  • $\begingroup$ Johan, it is absolutely fantastic having both you and @MichaelGrant helping out on scicomp. Do you mind ensuring that you disclosure your association with YALMIP in your answers? $\endgroup$
    – Aron Ahmadia
    Mar 20, 2013 at 18:05
  • $\begingroup$ Absolutely, sure $\endgroup$
    – Johan Löfberg
    Mar 20, 2013 at 18:33
  • $\begingroup$ This is becoming "meta" quickly, but if there were a way to add a line to our avatar disclosing our software ties... (i.e., "Michael Grant \\ CVX")... $\endgroup$
    – Michael Grant
    Mar 20, 2013 at 18:34
  • 1
    $\begingroup$ @MichaelGrant: That's a good meta question. Stack Exchange employees look at meta.scicomp.stackexchange.com and answer some of those questions. It's probably appropriate to apply the "feature-request" tag to such a question. $\endgroup$
    – Geoff Oxberry
    Mar 21, 2013 at 3:38
1
$\begingroup$

It's not clear to me what you're asking for here. Are you looking for a derivation of the semidefinite programming (SDP) formulation of the problem? Are you looking for an explanation of what SDP is? Are you looking for information on how to solve the SDP once you've got it?

MATLAB doesn't have any built-in functions for solving semidefinite programming problems. There are a number of open source packages that you can install that will add this capability to MATLAB. Look at SeDuMi, SDPT3, and CSDP. A quick google search will lead you to each of these packages.

Improve this answer
$\endgroup$
3
  • $\begingroup$ I'm still newly exposed to this linear programming stuff. I've updated the question, hopefully that helps. Basically my overarching quest is to find a method of estimating the maximum of a trig. polynomial and its location with a method that is faster than using MATLAB roots on the derivative of the poly. $\endgroup$
    – geometrikal
    Mar 20, 2013 at 0:16
  • $\begingroup$ I'm hoping to find a method that allows me to sacrifice accuracy for speed. $\endgroup$
    – geometrikal
    Mar 20, 2013 at 0:16
  • $\begingroup$ Wait, is it even possible to handle trigonometric polynomials with MATLAB's roots command? $\endgroup$
    – Michael Grant
    Mar 20, 2013 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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