%% 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 and updated in June 2026. CVX [2] is required to execute the optimization programs found here.

clear all; clc;

%
% The experiments below illustrate the optimal procedures uncovered in [1] for the estimation of a quadratic form \gamma associated with 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
gwce_min = cvx_optval;
% determination of the lower bound
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. For details about this computation, see the paragraph before Section 5.

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., a \not= 0 in
% \gamma(f) = \langle A f, f \rangle + 2 \langle a, f \rangle + \alpha.

clear all; clc;
% parameters of the problem
N = 50;
A_aux = rand(N,N);
A = (A_aux+A_aux')/2;
%a = rand(N,1)/sqrt(N);
a = rand(N,1);
alp = rand;
T = eye(N)/sqrt(N);
m = 30;
L = randn(m,N);

%
% For a fixed observation vector y \in \mathbb{R}^m, the minimal local worst-case error, aka Chebyshev radius, along with the Chebyshev center, are computed by solving two semidefinite programs (see Theorem 1). They are integrated in the function {\rm cheb} created at the end of this file.

f0_aux = randn(N,1);
f0 = 0.9*f0_aux/norm(T*f0_aux);
y = L*f0;
[cheb_center, cheb_radius] = cheb(A,a,alp,T,L,y);
fprintf(['The minimal local worst-case error is %.3f,\n'...
    'while the Chebyshev center is %.3f.'], cheb_radius, cheb_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, for any admissible y, y+h, y-h,
% q(y,h) := \delta^{\rm cheb}(y+h) + \delta^{\rm cheb}(y-h) - 2 \, \delta^{\rm cheb}(y)
% =
% 2 \langle Bh,h \rangle,
% so that the quantity q(y,h) would be independent of the choice of y. The experiment below refutes this independence by computing q(y_0,h) and q(y_1,h) for two different y_0 , y_1.

% define y0=L*f0,y1=L*f1 with f0,f1 strictly inside the model set
% define perturbations y0p=L*(f0+h),y0m=L*(f0-h),y1p=L*(f1+h),y1m=L*(f1-h)
% with h small enough for f0+h,f0-h,f1+h,f1-h to all belong to the model set
f0_aux = randn(N,1);
f0 = 0.5*f0_aux/norm(T*f0_aux);
y0 = L*f0;
f1_aux = rand(N,1);
f1 = 0.5*f1_aux/norm(T*f1_aux);
y1 = L*f1;
h_aux = randn(N,1);
h = 0.4*h_aux/norm(T*h_aux);
y0p = L*(f0+h);
y0m = L*(f0-h);
y1p = L*(f1+h);
y1m = L*(f1-h);
% compute the quantity q for y0
q0 = cheb(A,a,alp,T,L,y0p) + cheb(A,a,alp,T,L,y0m) - 2*cheb(A,a,alp,T,L,y0);
% compute the quantity q for y1
q1 = cheb(A,a,alp,T,L,y1p) + cheb(A,a,alp,T,L,y1m) - 2*cheb(A,a,alp,T,L,y1);
fprintf(['The above quantity q computed for y0 equals %.3f,\n'...
    'which differs from the one for y1, namely   %.3f.'], q0, q1)

%
% To conclude, it is also verified that the above claim (the locally optimal estimation functional is not a quadratic functional) is valid even when the quantity of interest is a quadratic form rather than a quadratic functional, i.e., when a=0.

a = zeros(N,1);
q0 = cheb(A,a,alp,T,L,y0p) + cheb(A,a,alp,T,L,y0m) - 2*cheb(A,a,alp,T,L,y0);
q1 = cheb(A,a,alp,T,L,y1p) + cheb(A,a,alp,T,L,y1m) - 2*cheb(A,a,alp,T,L,y1);
fprintf(['The above quantity q computed for y0 equals %.3f,\n'...
    'which differs from the one for y1, namely   %.3f.'], q0, q1)

%
%% 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.
%
%% Function computing the locally optimal recovery map
% The function {\rm cheb} defined below computes the Chebyshev center and radius 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 .

function [cheb_center, cheb_radius] = cheb(A,a,alp,T,L,y)
% auxiliary quantities
N = length(a);
invT = inv(T);
P = eye(N) - L'*inv(L*L')*L;
fy = (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*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;
cheb_radius = (sig-iot)/2;
cheb_center = (sig+iot)/2;
end