qpapproximation.m

%% Generate data
n = 10;
Q = randn(n);Q = Q*Q';
R = randn(n);R = R*R';
S = chol(Q);
T = chol(R);

%% Basic factorized model
x = sdpvar(n,1);
y = sdpvar(n,1);
e = sdpvar(n,1);
f = sdpvar(n,1);
Model = [-1 <= [x y] <= 1, sum(x) + sum(y) == 1, e == S*x, f == T*y];

%% Build approximation
[~,Le,Ue] = boundingbox(Model,[],e);
[~,Lf,Uf] = boundingbox(Model,[],f);
N = 100;
E = repmat(Le,1,N) + repmat(linspace(0,1,N),n,1).*repmat(Ue-Le,1,N);
F = repmat(Lf,1,N) + repmat(linspace(0,1,N),n,1).*repmat(Uf-Lf,1,N);
f1 = interp1(E,E.^2,e,'lp');
f2 = interp1(F,F.^2,f,'lp');

%% See what we are doing
z = -1:0.25:1;
mesh(z,z,z.^2-(z').^2)

%% Solve!
optimize(Model,sum(f1)-sum(f2))

%% Alternative with convex QP part
Model = [-1 <= [x y] <= 1, sum(x) + sum(y) == 1, f == T*y];
optimize(Model,x'*Q*x-sum(f2))


%% Example generic case
n = 10;
x = sdpvar(n,1);
y = sdpvar(n,1);
Q = randn(n);Q = Q*Q';
R = randn(n);R = R*R';
p = x'*Q*x - y'*R*y;

%% Assume we did not know structure
[H,c,b,z] = quaddecomp(p);
[V,D] = eig(full(H));
pos = find(diag(D)>0);
neg = find(diag(D)<0);
S = D(pos,pos)^.5*V(:,pos)';
T = (-D(neg,neg))^.5*V(:,neg)';

e = sdpvar(size(S,1),1);
f = sdpvar(size(T,1),1);
Model = [-1 <= [x y] <= 1, sum(x) + sum(y) == 1, e == S*z, f == T*z];
[~,Le,Ue] = boundingbox(Model,[],e);
[~,Lf,Uf] = boundingbox(Model,[],f);
N = 100;
E = repmat(Le,1,N) + repmat(linspace(0,1,N),n,1).*repmat(Ue-Le,1,N);
F = repmat(Lf,1,N) + repmat(linspace(0,1,N),n,1).*repmat(Uf-Lf,1,N);
f1 = interp1(E,E.^2,e,'lp');
f2 = interp1(F,F.^2,f,'lp');
optimize(Model,sum(f1)-sum(f2) + c'*z + b)

%% With convex part kept
Model = [-1 <= [x y] <= 1, sum(x) + sum(y) == 1, f == T*z];
optimize(Model,z'*S'*S*z - sum(f2) + c'*z + b)