%% Reproducible file accompanying the paper
%% WORST-CASE-OPTIMAL ESTIMATION OF A QUADRATIC FORM BY QUADRATIC FUNCTIONALS OF ITS LINEAR OBSERVATIONS
% by S. Foucart.
%
% This reproducible file was written by S. Foucart in March 2026. CVX [2] is required to execute the optimization programs found here. Note that some random choices of parameters may result in semidefinite programs that do not solve.

clear all; clc;

%
% The experiments below illustrate the optimal procedures uncovered in [1] for the estimation of a quadratic form \gamma associated to a self-adjoint operator A on a Hilbert space H, namely
% \gamma(f) = \langle Af, f \rangle ,
% given partial information on the unknown f \in H, precisely some a priori information in the form
% \|T f \| \le 1, \qquad T:H \to H \mbox{ linear operator},
% and some a posteriori information in the form
% y = \Lambda f,
% \qquad
% \Lambda: H \to \mathbb{R}^m \mbox{ linear map}.
%
%% Setting of the problem
% The problem parameters are introduced here: the dimension N of the underlying Hilbert space (H = \ell_2^N(\mathbb{R}) here), the quadratic form \gamma (via a symmetric matric A), the operator T, and the observation map \Lambda (denoted by L below).

N = 50;
A_aux = rand(N,N);
A = (A_aux+A_aux')/2;
T = randn(N,N)/sqrt(N);
m = 30;
L = randn(m,N);

%
%% Globally Optimal Estimation Functionals: Accurate Scenario
% One considers first the global setting, where the performance of an estimation functional \delta: \mathbb{R}^m \to \mathbb{R} is assessed, in the absence of observation errors, through
% {\rm gwce}(\delta) = \sup \{ |\delta(f) - \Delta(\Lambda f)| : \|Tf \| \le 1\}.
% One starts by determining the minimal global worst-case error over all estimation functionals, computed according to Theorem 6. Recall that this relies on the fact that quadratric estimation functionals are optimal, as established in Theorem 2. This fact is also corroborated empirically by verifying the numerical agreement with a general lower bound, computed according to Proposition 8.

% minimization over quadratic functionals
cvx_begin quiet
variable cp nonnegative
variable cpp nonnegative
variable B(m,m) symmetric
minimize .5*(cp+cpp)
subject to
cp*T'*T  - A + L'*B*L == semidefinite(N);
cpp*T'*T + A - L'*B*L == semidefinite(N);
cvx_end
% determination of the lower bound
gwce_min = cvx_optval;
cvx_begin quiet
variable c1 nonnegative
variable c2 nonnegative
variable d
minimize .5*(c1+c2)
subject to
[c1*T'*T - A + d*L'*L, -d*L'*L;...
    -d*L'*L, c2*T'*T + A + d*L'*L]...
    == semidefinite(2*N);
cvx_end
lb = cvx_optval;
fprintf(['The minimal global worst-case error'...
    ' and the general lower bound agree.\n' ...
    'They are (up to numerical errors) %.3f and %.3f, respectively.'],...
    gwce_min,lb)

%
% As a next step, one inquires about the global worst-case error of the "plug-in" estimation functional, i.e., of \delta^{\rm plug} := \gamma \circ \Delta^{\rm opti}  (plus a scalar), where \Delta^{\rm opti}: \mathbb{R}^m \to H    is the locally (hence globally) optimal map for the recovery of f in full.

invT = inv(T);
Del_opti = (invT*invT')*L'*inv(L*(invT*invT')*L');
cvx_begin quiet
variable cp nonnegative
variable cpp nonnegative
minimize .5*(cp+cpp)
subject to
cp*T'*T  - A + L'*Del_opti'*A*Del_opti*L == semidefinite(N);
cpp*T'*T + A - L'*Del_opti'*A*Del_opti*L == semidefinite(N);
cvx_end
gwce_plug = cvx_optval;
fprintf(['The global worst-case error of the plug-in functional equals %.3f, \n'...
    'compared with the minimal global worst-case error, which was %.3f.'],gwce_plug,gwce_min)

%
% The above experiment, conducted with a "random" operator T, can lead to the belief that the "plug-in" estimation functional is optimal in general. However, this is not the case, as revealed by the experiment below conducted with T = {\rm Id}_H.

% minimization over quadratic functionals
cvx_begin quiet
variable cp nonnegative
variable cpp nonnegative
variable B(m,m) symmetric
minimize .5*(cp+cpp)
subject to
cp*eye(N)  - A + L'*B*L == semidefinite(N);
cpp*eye(N) + A - L'*B*L == semidefinite(N);
cvx_end
gwce_min_id = cvx_optval;
% determination of the global worst-case error of the plug-in 
Del_opti_id = L'*inv(L*L');
cvx_begin quiet
variable cp nonnegative
variable cpp nonnegative
minimize .5*(cp+cpp)
subject to
cp*eye(N)  - A + L'*Del_opti_id'*A*Del_opti_id*L == semidefinite(N);
cpp*eye(N) + A - L'*Del_opti_id'*A*Del_opti_id*L == semidefinite(N);
cvx_end
gwce_plug_id = cvx_optval;
fprintf(['The global worst-case error of the plug-in functional equals %.3f, \n'...
    'while the minimal global worst-case error is now equal to %.3f.'],...
    gwce_plug_id,gwce_min_id)

%
%% Globally Optimal Estimation Functionals: Inaccurate Scenario
% In the presence of deterministic observation errors e \in \mathbb{R}^m with \|e\|_2 \le \eta, the global worst-case error of an estimation functional \delta: \mathbb{R}^m \to \mathbb{R} is adjusted to become
% {\rm gwce}(\delta) = \sup \{ |\delta(f) - \Delta(\Lambda f + e)| : \|Tf \| \le 1, \|e\|_2 \le 1\}.
% Theorem 7 guarantees that the minimum over all estimation functionals is again achieved for a quadratic functional and that it can still be computed by solving a semidefinite program. Note that it has to be larger than the minimal global worst-case error in the accurate setting.

U = eye(m);
eta = 0.01;
cvx_begin quiet
variable cp1 nonnegative
variable cp2 nonnegative
variable cpp1 nonnegative
variable cpp2 nonnegative
variable B(m,m) symmetric
minimize .5*(cp1+cpp1+eta^2*(cp2+cpp2))
subject to
[cp1*T'*T - A + L'*B*L, L'*B;...
    B*L, cp2*U'*U + B]...
    == semidefinite(N+m);
[cpp1*T'*T + A - L'*B*L, -L'*B;...
    -B*L, cpp2*U'*U - B]...
    == semidefinite(N+m);
cvx_end
gwce_min_inac = cvx_optval;
fprintf(['The minimal global worst-case error in the inaccurate scenario is %.3f, \n'...
    'compared with the minimal one in the accurate scenario, which was %.3f.'],gwce_min_inac,gwce_min)

%
%% Locally Optimal Estimation Functionals
% Given some fixed y \in \mathbb{R}^m obtained as y = \Lambda f_0 for some f_0 \in H satisfying \|T f_0 \| \le 1, the objective is now to find the Chebyshev center of the set \{ \gamma(f): f \in H \mbox{ satisfies } \|Tf \| \le 1, \Lambda f = y\}, i.e.,
% g^{\rm cheb} = \underset{g \in \mathbb{R}}{\rm argmin}
% \sup_{f \in H} \{ |\gamma(f) - g|: \|Tf\| \le 1, \Lambda f = y \}.
% The considerations are extended to the estimation of quadratic functionals instead of merely quadratic forms, i.e.,
% \gamma(f) = \langle A f, f \rangle + 2 \langle a, f \rangle + \alpha.
% The Chebyshev center is computed according to Theorem 1, whose statement involves the orthogonal projector P onto \ker(\Lambda) and the element f_y \in H minimizing \|Tf\| subject to \Lambda f = y .

a = rand(N,1)/N;
alp = rand;
T = eye(N)/sqrt(N);
f0_aux = randn(N,1);
f0 = 0.9*f0_aux/norm(T*f0_aux);
y = L*f0_aux;
P = eye(N) - L'*inv(L*L')*L;
fy = (invT*invT')*L'*inv(L*(invT*invT')*L')*y;

%
% The minimal local worst-case error is computed, along with the Chebyshev center, by solving two semidefinite programs.

% solve the first SDP
cvx_begin quiet
variable c nonnegative
variable d
minimize d
subject to
[P*(c*(T'*T) - A)*P,...
    -P*(A*fy+a); -(A*fy+a)'*P,...
    d - c*(1-norm(T*fy)^2) - fy'*A*fy - 2*a'*fy - alp]...
    == semidefinite(N+1);
cvx_end
sig = cvx_optval;
% solve the second SDP
cvx_begin quiet
variable c nonnegative
variable d
maximize d
subject to
[P*(c*(T'*T) + A)*P,...
    P*(A*fy+a); (A*fy+a)'*P,...
    -d - c*(1-norm(T*fy)^2) + fy'*A*fy + 2*a'*fy + alp]...
    == semidefinite(N+1);
cvx_end
iot = cvx_optval;
lwce = (sig-iot)/2;
center = (sig+iot)/2;
fprintf(['The minimal local worst-case error is %.3f,\n'...
    'while the Chebyshev center is  %.3f.'], lwce, center)

%
% It can be verified empirically that the locally optimal estimation functional \delta^{\rm cheb}: y \in \mathbb{R}^m \mapsto g^{\rm cheb} \in \mathbb{R} is not quadratic. Indeed, if it had the form
% \delta^{\rm cheb}(y) = \langle By, y \rangle + 2 \langle b,y \rangle + \beta,
% \qquad y \in \mathbb{R}^m,
% then one would have \delta^{\rm cheb}(y) - \delta^{\rm cheb}(-y) = 4 \langle b, y \rangle for all admissible y, and in turn
% \delta^{\rm cheb}(y) - \delta^{\rm cheb}(-y) - 2 (\delta^{\rm cheb}(y/2) - \delta^{\rm cheb}(-y/2)) = 0.
% Computing these four quantities shows that the above expression is not equal to zero.
%

% for the first delta, use the previous computation
del_y = center;
% for the second delta, replace y by the opposite of y
fyopp = -(invT*invT')*L'*inv(L*(invT*invT')*L')*y;
% solve the first SDP
cvx_begin quiet
variable c nonnegative
variable d
minimize d
subject to
[P*(c*(T'*T) - A)*P,...
    -P*(A*fyopp+a); -(A*fyopp+a)'*P,...
    d - c*(1-norm(T*fyopp)^2) - fyopp'*A*fyopp - 2*a'*fyopp - alp]...
    == semidefinite(N+1);
cvx_end
sig = cvx_optval;
% solve the second SDP
cvx_begin quiet
variable c nonnegative
variable d
maximize d
subject to
[P*(c*(T'*T) + A)*P,...
    P*(A*fyopp+a); (A*fyopp+a)'*P,...
    -d - c*(1-norm(T*fyopp)^2) + fyopp'*A*fyopp + 2*a'*fyopp + alp]...
    == semidefinite(N+1);
cvx_end
iot = cvx_optval;
del_yopp = (sig+iot)/2;
% for the third delta, replace y by the half of y
fyhalf = (invT*invT')*L'*inv(L*(invT*invT')*L')*y/2;
% solve the first SDP
cvx_begin quiet
variable c nonnegative
variable d
minimize d
subject to
[P*(c*(T'*T) - A)*P,...
    -P*(A*fyhalf+a); -(A*fyhalf+a)'*P,...
    d - c*(1-norm(T*fyhalf)^2) - fyhalf'*A*fyhalf - 2*a'*fyhalf - alp]...
    == semidefinite(N+1);
cvx_end
sig = cvx_optval;
% solve the second SDP
cvx_begin quiet
variable c nonnegative
variable d
maximize d
subject to
[P*(c*(T'*T) + A)*P,...
    P*(A*fyhalf+a); (A*fyhalf+a)'*P,...
    -d - c*(1-norm(T*fyhalf)^2) + fyhalf'*A*fyhalf + 2*a'*fyhalf + alp]...
    == semidefinite(N+1);
cvx_end
iot = cvx_optval;
del_yhalf = (sig+iot)/2;
% for the fourth delta, replace y by the half of the opposite of y
fyopphalf = -(invT*invT')*L'*inv(L*(invT*invT')*L')*y/2;
% solve the first SDP
cvx_begin quiet
variable c nonnegative
variable d
minimize d
subject to
[P*(c*(T'*T) - A)*P,...
    -P*(A*fyopphalf+a); -(A*fyopphalf+a)'*P,...
    d - c*(1-norm(T*fyopphalf)^2) - fyopphalf'*A*fyopphalf - 2*a'*fyopphalf - alp]...
    == semidefinite(N+1);
cvx_end
sig = cvx_optval;
% solve the second SDP
cvx_begin quiet
variable c nonnegative
variable d
maximize d
subject to
[P*(c*(T'*T) + A)*P,...
    P*(A*fyopphalf+a); (A*fyopphalf+a)'*P,...
    -d - c*(1-norm(T*fyopphalf)^2) + fyopphalf'*A*fyopphalf + 2*a'*fyopphalf + alp]...
    == semidefinite(N+1);
cvx_end
iot = cvx_optval;
del_yopphalf = (sig+iot)/2;
expr = del_y - del_yopp - 2*(del_yhalf - del_yopphalf);
fprintf(['If the locally optimal estimation functional was quadratic,\n'...
    'the following value would be zero: %.3f.'], expr)

%
%% References
% S. Foucart, "Worst-case-optimal estimation of a quadratic form by quadratic functionals of its linear observations", Preprint.
% M. Grant and S. Boyd, "CVX: Matlab software for disciplined convex programming", version 2.1. http://cvxr.com/cvx, March 2014.