Minimizing the values in an underdetermined Matrix - minimize with constraints

I'm new to MATLAB, and there is probably a very simple way to do this. I have a linear set of equations Ax=b. I have matrix A and matrix b, and I need to solve for the smallest values in matrix x that satisfy the equation. The solution is underdetermined. Is there a way to do this? Essentially I need to minimize the components in the matrix x subject to constraints. Mathematica has the NMinimize function which looks perfect, except I cannot find an equivalent in MATLAB. fmincon does not seem to work for matrix entries.

Thanks

[Richard Brown]

Out of interest, are you doing compressed sensing? l1-magic is a collection of MATLAB routines designed for such problems: http://users.ece.gatech.edu/~justin/l1magic/

Answer 1

[Richard Brown]

edit: fixed mistake in call to linprog

To perform L1 minimisation, you'll need access to an LP solver. If you have the optimization toolbox, linprog is your friend.

The easiest way to do it is as follows:

Define a vector t of the same length as x such that -t <= x <= t

The linear program is then

minimise    t(1) + ... + t(n)
subject to  A*x = b
            -x - t <= 0
             x - t <= 0

To solve it in MATLAB, I'll assume you have an m x n matrix A, and m x 1 vector b

[m, n] = size(A); 
f = [zeros(n, 1); ones(n, 1)];

Ai = [-eye(n), -eye(n); eye(n), -eye(n)];
bi = zeros(2*n, 1);

x = linprog(f, Ai, bi, [A, zeros(m, n)], b);
x = x(1:n);

[Joseph]

Getting an error; "number of columns in Aeq must be the same as the length of f." Not sure if this is because it is an underdetermined system. I am trying this with A = 5x10 random matrix, b = 5x1 matrix. I have tried changing a few things, no luck yet.

[Richard Brown]

Sorry, you need to replace the A in the call to linprog, by [A, zeros(m, n)]

[Joseph]

Works Great! Thanks.

[Joseph]

Is there a way to solve this for a complex case? When A and b are complex? Thanks, Joseph

[Giorgio]

Thank you for this helpful code! How should I modify the matrix Ai in order not to have negative results (x>=0)?

[Torsten]

[m, n] = size(A); 
f = [zeros(n, 1); ones(n, 1)];

Ai = [-eye(n), -eye(n); eye(n), -eye(n)];
bi = zeros(2*n, 1);

lb = zeros(2*n, 1);
ub = Inf(2*n, 1);

x = linprog(f, Ai, bi, [A, zeros(m, n)], b, lb, ub);
x = x(1:n);

Best wishes

Torsten.

Answer 2

[Richard Brown]

Assuming you want a minimal Euclidean norm solution, there's no built-in function, but I can offer you a two-line solution:

[Q, R] = qr(A', 0);
x = Q * (R' \ b);

If you want a minimal 1 or infinity norm solution, then you can cast the problem as a linear program.

[Joseph]

I actually need the Taxicab or L1 norm. Can you explain how your above code works?

[Teja Muppirala]

That's a neat trick, Richard! Instead of explicitly calling QR though, I think the following variants may be slighty faster:

x = ((A*A')\(A))'*b;

Or

x = (A'/(A*A'))*b;

I would disagree that there are no built-in solutions though. There are many general purpose optimizers available for MATLAB, which while may be overkill for this specific problem, are certainly capable of solving it (and other more difficult problems) nearly as easily. For example, in this case, even if one did not come up with the simple analytic solution that you presented, you could have still solved it like either of these:

x = quadprog(eye(length(A)),zeros(1,length(A)),[],[],A,b)

x = fmincon(@(x) norm(x), zeros(1,length(A)),[],[],A,b)

[Richard Brown]

@Joseph (btw to respond to an answer, use comments underneath the question rather than a new answer)

To understand the code above, consider the full QR factorization of A'

[Qfull, Rfull] = qr(A')

The orthogonal matrix Qfull can be partitioned as [Q Qnull] and Rfull as [R; 0]. Define y = Qfull' * x and partition it as y = [u; v]. Then, Ax = b becomes

Rfull' * Qfull' * x = b
Rfull' * y = b
R' * u = b

So u = R' \ b and v is free to choose.

Now norm(x, 2)^2 = norm(y, 2)^2 = norm(u)^2 + norm(v)^2. Since v is free to choose, choose v = 0 to minimise the norm. Then

|x = Qfull * y = Q *u + Qnull * v 
   = Q * (R' \ b)|.

So essentially you're minimising norm(x) over the null space of A. Finally, you don't need to calculate the full QR factorization, you can get Q and R by

[Q, R] = qr(A', 0);

[Richard Brown]

@Teja, you're right, it's definitely faster, because MATLAB can do a Cholesky factorisation instead of a QR. By the way, it's much faster if you rewrite it as

x = A' * ((A*A')\ b);

(note to those interested - this is the solution you get if you enforce x to be in the range of A', i.e. no component in the nullspace of A. In this case, set x = A' * y. Then

A * A' y = b 
=> y = (A * A') \ b
=> x = A' * ((A * A') \ b)

)

BUT your method is much less stable - the matrix A*A' can be very badly conditioned, so it's not recommended!

And yes, I know you can use the optimization toolbox routines to solve it. I meant no specific builtin. This solution is what I think A \ b should return in this case.

[Joseph]

x = fmincon(@(x) norm(x,1), zeros(1,length(A)),[],[],A,b)

The above worked for me. Thanks Teja. I actually tried the above but without the @x - What does the @x do?

[Richard Brown]

Glad that worked for you. When I tried it on some random examples (even small ones), your call to fmincon with default parameters didn't converge. This isn't surprising, because the 1-norm objective isn't differentiable, and fmincon is designed for smooth problems.

I would still recommend using an LP based solution if you're interested in performance - I'd expect you'd get 10x or better speed performance from linprog for your problem