MATLAB: Reformulate a Constrained Linear Least Square Problem

constrained optimizationslinear optimizationslsqlinMATLABOptimization Toolbox

I have a problem of the form
Normally you solve it like this:
u = (0:0.1:10)';
v = sin(2*pi*1/30*(u-0.5));
plot(u,v);
xlabel('u');
ylabel('v');
C = u.^[0 1 2 3 4 5 6 7];
d = v;
u2 = (0:0.01:11)';
A = -u2.^[0 1 2 3 4 5 6 7];
b = zeros(size(A,1),1);
% normally just use lsqlin to solve this
x = lsqlin(C,d,A,b);
hold all
plot(u,C*x)
However, my problem is actually quite large and ill-conditioned because of the high degree of polynomial that I am fitting. The interior-point solver uses C'*C to solve this problem. That is a problem because it squares the condition number of C. The interior-point algorithm struggles to converge while the older, now deprecated Active-set algorithm works well. The Active-set algorithm wasn't large scale. It was medium scale. It worked great. However, I don't have that available. What are some ways I can stabilize and solve this problem?
Some ideas:
  • Maybe I can use the null space of A?
  • Maybe I can reformulate the problem to fmincon with proper settings.
  • I know about orthogonal transformation for a single dimension polynomial. However, I work with multidimensional polynomials and even though orthogonal decomposition and formulations exist, I don't know them or how to use them. If you propose this, please refer to how I can use these and where I can learn about them and understand them.

Best Answer

You can try a trick from Lawson & Hanson and re-formulate to a minimal distance problem and solve via lsqnonneg. lsqnonneg is an active-set method, so if those work for you, then it might as well.
The minimum distance problem looks like:
   minimize ||x||^2
      x
     s.t.   Abar*x <= bbar
Code:
% Transform into minimal distance
[Q,R] = qr(C,0);
dbar = Q'*d;
Abar = A/R;
bbar = b - ARinv*dbar;
% Get min-distance into lsqnonneg form 
n = size(Abar,2);
E = [Abar'; 
     bbar'];
f = [zeros(n,1); -1];
[u,~,residual] = lsqnonneg(E,f);
xbar = -residual(1:n)/residual(end);
% Map back
x = R\(xbar+dbar);
Related Question