|
Can YALMIP help solve this problem? 1 Year, 9 Months ago
|
Karma: 0
|
|
Hello all,
I hope you can advise me if YALMIP can help with the following problem.
min I1(N,M) (complex) , r(N,M) (real)
objective = ||data1 - F1*I1*F2||2 + ||data2 - F1*(I1.*r)*F2||2 + lambda*( ||Dx*r||2 + ||Dy*r||2 )
In context, I have formulated a problem in medical image reconstruction as above. I1 is a complex image and r is a 'smooth' image. The first two terms are data consistency terms, and the last two are spatial finite different operators.
The problem as I see it is the second term, where I1 and r interact, making the problem nonlinear in the optimisation variables. Unfortunately, my background in optimisation is not sufficient to formulate this problem in a tractable way!
Thanks for you help, it's very much appreciated.
|
|
|
|
|
|
|
Re: Can YALMIP help solve this problem? 1 Year, 9 Months ago
|
Karma: 32
|
Basically no (well, it can model it, YALMIP is a language and not primarily a solver)
Due to the products between l1 and r you have a pretty nasty nonconvex problem.
Since you are talking about an image too, I guess the problem is large, so I would guess you best choice is a dedicated local nonlinear programming approach
a trivial iterative SOCP scheme is simple to implement though
| Code: |
l1 = sdpvar(n,m,'full','complex');
r = sdpvar(n,m,'full');
% First guess
assign(r,zeros(n,m));
while 1
objective = norm(data1 - F1*I1*F2,2)+norm(data2 - F1*(I1.*double(r))*F22,2)+lambda*(norm(Dx*double(r))+norm(Dy*double(r)));
solvesdp([],objective);
objective = norm(data1 - F1*double(I1)*F2,2)+norm(data2 - F1*(double(I1).r)*F22,2)+lambda*(norm(Dx*r)+norm(Dy*r));
solvesdp([],objective);
end
|
|
|
lofberg
Platinum Boarder
Posts: 2280
|
|
Last Edit: 2011/08/16 06:40 By lofberg.
|
|
|
Re: Can YALMIP help solve this problem? 1 Year, 9 Months ago
|
Karma: 0
|
|
Thank you profuseley for the prompt reply, Lofberg!
The problem isn't huge - N = 32,M = 20, but I guess it depends on what you're used to. Do you have any pointers for the local nonlinear programming methods I could look at?
As for the code you wrote, I am trying this at the moment using CVX - however, I think there are no guarentees on convergence. I'm also trying NL conjugate gradients.
Perhaps a genetic-type algorithm could work?
Again, thanks for the help!
|
|
|
|
|
|
|
Re: Can YALMIP help solve this problem? 1 Year, 9 Months ago
|
Karma: 32
|
|
cvx vs yalmip here doesn't matter. In the end, you will have to come with some scheme which solves SOCPs repeatedly. Just a matter of different modeling of your algorithm.
I have no direct reference, but it looks like a problem (or variant of) which has been addressed in various ways in the image processing community. I am not expert in that field though
A standard vanilla gradient based approach might have problems, since you're using non-smooth objectives.
Genetic-type seems like the wrong approach, problem looks to structured.
Are you married to the use of sum of norms instead of sum of squared norms? Sum of squared norms would probably be a lot easier (smooth etc), and you go for an iterative scheme, you an iteration of least-squares problem, i.e. trivial.
|
|
lofberg
Platinum Boarder
Posts: 2280
|
|
|
|
|
Re: Can YALMIP help solve this problem? 1 Year, 9 Months ago
|
Karma: 0
|
|
Sorry for the confusion, but I do currently use squared norms - I simply forgot to mention in the OP. This makes the function smooth, as you say, but the nonlinearity still leads to problems, as far as I'm aware. I've seen solutions fall into local minima, and erratic solution behavior for small changes in smoothness regularization parameter. Hence the desire to use CVX or Yalmip, 'trust-worthy' toolboxes which would hopefully provide reliable solutions.
Hope these developments give you some more insights into the problem!
|
|
|
|
Last Edit: 2011/08/16 08:34 By Antonius1.
|
|
|
Re: Can YALMIP help solve this problem? 1 Year, 9 Months ago
|
Karma: 32
|
Well, if the norms you are using are frobenious norms on the complex-valued matrices, the problem boils down to a multilinear least-squares problem which is a well-known hard problem.
Hence, when you want to solve this problem for real, you shouldn't use a modelling layer such as cvx or YALMIP, as you can code blazingly fast algorithms for this problem (simple derivatives, based on least squares etc).
this would be the YALMIP version of the mulltilinear least-squares, i.e. each iteration is a trivial QP, and not an SOCP which you get using cvx, or YALMIP with norms.
note that 99% of computations are spent in the modelling, the QP is solved in 0 seconds.
| Code: |
n = 5;
m = 6;
k = 7;
I1 = sdpvar(n,m,'full','complex');
r = sdpvar(n,m,'full');
data1 = randn(k,k)+randn(k,k)*sqrt(-1);
data2 = randn(k,k);
F1 = randn(k,n);
F2 = randn(m,k);
F22 = randn(m,k);
Dx = eye(n);
Dy = eye(n);
lambda = 1;
assign(r,zeros(n,m));
while 1
objective = norm(data1 - F1*I1*F2,'fro')^2+norm(data2 - F1*(I1.*double(r))*F22,'fro')^2+lambda*(norm(Dx*double(r),'fro')^2+norm(Dy*double(r),'fro')^2);
solvesdp([],objective);
double(objective)
objective = norm(data1 - F1*double(I1)*F2,'fro')^2+norm(data2 - F1*(double(I1).*r)*F22,'fro')^2+lambda*(norm(Dx*r,'fro')^2+norm(Dy*r,'fro')^2);
solvesdp([],objective);
double(objective)
end
|
|
|
lofberg
Platinum Boarder
Posts: 2280
|
|
|
|
|
