%% Reproducible file accompanying the paper
%% OPTIMAL ALGORITHMS FOR NONLINEAR ESTIMATION WITH CONVEX MODELS
% by S. Foucart.
% This reproducible was written by S. Foucart in December 2025.
% CVX [2] is required to execute this reproducible.
%
%% Estimation of the maximum of linear functionals with polytopal models.

clear all; clc;

%
% As descrived in Subsection 3.2 of [1], in the space F = \ell_\infty^N, one considers the polytopal model set
% \mathcal{K}^{\rm poly} = \{ f \in \ell_\infty^N: \langle a_\ell, f \rangle \le 1 \mbox{ for } \ell \in L \}.
% The linear observations \Lambda(f) \in \mathbb{R}^m made on an unknown f \in \mathcal{K}^{\rm poly} are, entrywise,
% \lambda_k(f) = \langle u_k, f \rangle,
% \qquad k=1,\ldots,m.
% The objective is to estimate the sup-linear functional \gamma defined by
% \gamma(f) = \sup_{i \in I} \gamma_i(f),
% \qquad \mbox{with} \quad
% \gamma_i(f) = \langle w_i, f \rangle,
% in a way that minimizes the worst-case error
% e(\delta) = \sup_{f \in \mathcal{K}^{\rm poly}} | \gamma(f) - \delta(\Lambda(f))|
% over all possible functionals \delta: \mathbb{R}^m \to \mathbb{R}. Here is the problem setting.

% introduce the model set 
% (may need to increase L_size to make the model set bounded)
N = 10;  
L_size = 30;
A = randn(N,L_size);
% introduce the observation process
m = 6;
U = randn(N,m);
% introduce the functional to be estimated
I_size = 8;
W = randn(N,I_size); 

%
% Theorem 1 reveals that an optimal estimation functional \delta^{\rm opti}: \mathbb{R}^m \to \mathbb{R} has a sup-affine form. It is determined by solving the following linear program stemming from Proposition 5.

cvx_begin quiet
variable e_opti
variable e1(I_size)
variable e2(I_size)
variable c(m,I_size)
variable sig(I_size,I_size) nonnegative
variable s(L_size,I_size) nonnegative
variable t(L_size,I_size) nonnegative
minimize e_opti
subject to 
sum(sig) == 1;
for i=1:I_size
    e1(i)+e2(i) <= 2*e_opti;
    sum(s(:,i)) <= e1(i);
    sum(t(:,i)) <= e2(i);
    A*s(:,i) == W(:,i) - U*c(:,i);
    A*t(:,i) == -W*sig(:,i) + U*c(:,i);
end
cvx_end
fprintf(['The minimal error over all estimation functionals '...
    'is equal to %.3f.'],e_opti)

%
% For the sake of verification, the optimal estimation functional \delta^{\rm opti}: \mathbb{R} \to \mathbb{R} is constructed and its error is computed via the linear program stemming from Proposition 4.

Z = c;
b = e_opti - e2; 
cvx_begin quiet
variable e_verif
variable sig(I_size,I_size) nonnegative
variable tau(I_size,I_size) nonnegative
variable s(L_size,I_size) nonnegative
variable t(L_size,I_size) nonnegative
minimize e_verif
subject to
sum(tau) == 1;
sum(sig) == 1;
for i=1:I_size
    sum(s(:,i)) <= e_verif + b'*sig(:,i);
    sum(t(:,i)) <= e_verif - b(i);
    A*s(:,i) == W(:,i) - U*Z*sig(:,i);
    A*t(:,i) == U*Z(:,i) - W*tau(:,i);
end
cvx_end
fprintf(['The verified error was computed to be %.3f, \n'...
    'indeed agreeing with the minimal error of %.3f.'],...
    e_verif,e_opti)

%
% Next, one constructs the optimal estimation map \Delta^{\rm opti}: \mathbb{R}^m \to \ell_\infty^N for the full approximation problem, which is an affine map whose components are computed by solving the following linear program.

cvx_begin quiet
variable e_full
variable c(m,N)
variable c0(N)
variable s(L_size,N) nonnegative
variable t(L_size,N) nonnegative
minimize e_full
subject to
sum(s) <= e_full + c0';
sum(t) <= e_full - c0';
A*s == eye(N) - U*c;
A*t == -eye(N) + U*c;
cvx_end

%
% One verifies that the plug-in functional for the estimation of \gamma, i.e., the natural candidate \delta^{\rm plug} := \gamma \circ \Delta^{\rm opti}: \mathbb{R}^m \to \mathbb{R}, leads to an error e(\delta^{\rm plug})---computed via the linear program stemming from  Proposition 4---larger than the minimal error e(\delta^{\rm opti}).

Z = c*W;
b = W'*c0;
cvx_begin quiet
variable e_plug
variable sig(I_size,I_size) nonnegative
variable tau(I_size,I_size) nonnegative
variable s(L_size,I_size) nonnegative
variable t(L_size,I_size) nonnegative
minimize e_plug
subject to
sum(sig) == 1;
sum(tau) == 1;
for i=1:I_size
    sum(s(:,i)) <= e_plug + b'*sig(:,i);
    sum(t(:,i)) <= e_plug - b(i);
    A*s(:,i) == W(:,i) - U*Z*sig(:,i);
    A*t(:,i) == U*Z(:,i) - W*tau(:,i);
end
cvx_end
fprintf(['The error of the plug-in functional is %.3f, \n'...
    'but the minimal error was smaller as %.3f.'],...
    e_plug,e_opti)

%
%% Estimation of the maximum of linear functionals with approximability models.

clear all; clc;

%
% As descrived in Subsection 3.3 of [1], in an N-dimensional Hilbert space F=H, one considers the approximability model set
% \mathcal{K}^{\rm appr} := \{ f \in H: \|P_{V^\perp}(f-g) \| \le \epsilon \}
% relative to a parameter \epsilon > 0, an n-dimensional linear subspace V \subseteq H, and an element g \in H \setminus V satisfying \|P_{V^\perp}g\| < \epsilon (so that 0 \in \mathcal{K}^{\rm appr}). Again, the linear observations \Lambda(f) \in \mathbb{R}^m made on an unknown f \in \mathcal{K}^{\rm appr} are, entrywise,
% \lambda_k(f) = \langle u_k, f \rangle,
% \qquad k=1,\ldots,m.
% The objective is once more to estimate the sup-linear functional \gamma defined by
% \gamma(f) = \sup_{i \in I} \gamma_i(f),
% \qquad \mbox{with} \quad
% \gamma_i(f) = \langle w_i, f \rangle,
% in a way that minimizes the worst-case error
% e(\delta) = \sup_{f \in \mathcal{K}^{\rm appr}} | \gamma(f) - \delta(\Lambda(f))|
% over all possible functionals \delta: \mathbb{R}^m \to \mathbb{R}. Here is the problem setting.

% introduce the model set 
N = 50; 
n = 20; 
V = rand(N,n);
P = eye(N) - V*inv(V'*V)*V';
epsilon = 0.7;
g_aux = rand(N,1);
g = (epsilon/2)*g_aux/norm(P*g_aux);
% introduce the observation process
m = 30;
U = randn(N,m);
U = U*inv(sqrtm(U'*U));      % to guarantee that U'*U=I
% introduce the functional to be estimated
I_size = 60;
W = randn(N,I_size); 

%
% Theorem 1 reveals that an optimal estimation functional \delta^{\rm opti}: \mathbb{R}^m \to \mathbb{R} has a sup-affine form. It is determined by solving the following second-order cone program stemming from Proposition 5.

cvx_begin quiet
variable e_opti
variable e1(I_size)
variable e2(I_size)
variable c(m,I_size)
variable sig(I_size,I_size) nonnegative
minimize e_opti
subject to 
sum(sig) == 1;
for i=1:I_size
    e1(i)+e2(i) <= 2*e_opti;
    V' * (W(:,i)-U*c(:,i)) == 0;
    V' * (-W*sig(:,i)+U*c(:,i)) == 0;
    epsilon*norm(W(:,i)-U*c(:,i)) <= ...
        e1(i) - g'*(W(:,i)-U*c(:,i));
    epsilon*norm(-W*sig(:,i)+U*c(:,i)) <= ...
        e2(i) - g'*(-W*sig(:,i)+U*c(:,i));
end
cvx_end
fprintf(['The minimal error over all estimation functionals '...
    'is equal to %.3f.'],e_opti)

%
% For the sake of verification, the optimal estimation functional \delta^{\rm opti}: \mathbb{R} \to \mathbb{R} is constructed and its error is computed via the second-order cone program stemming from Proposition 4.

Z = c;
b = e_opti - e2;
cvx_begin quiet
variable e_verif
variable sig(I_size,I_size) nonnegative
variable tau(I_size,I_size) nonnegative
minimize e_verif
subject to
sum(sig) == 1;
sum(tau) == 1;
for i=1:I_size
    V' * (W(:,i)-U*Z*sig(:,i)) == 0;
    V' * (U*Z(:,i)-W*tau(:,i)) == 0;
    epsilon*norm(W(:,i)-U*Z*sig(:,i)) <= ...
        e_verif + b'*sig(:,i) - g'*(W(:,i)-U*Z*sig(:,i)); 
    epsilon*norm(U*Z(:,i)-W*tau(:,i)) <= ...
        e_verif - b(i) - g'*(U*Z(:,i)-W*tau(:,i));
end
cvx_end
fprintf(['The verified error was computed to be %.3f, \n'...
    'indeed agreeing with the minimal error of %.3f.'],...
    e_verif,e_opti)

%
% Next, one constructs the optimal estimation map \Delta^{\rm opti}: \mathbb{R}^m \to H for the full approximation problem. It is an affine map which has the form
% \Delta^{\rm opti}(y)
% = g - \Delta^{\rm cheb}(\Lambda(g)) + \Delta^{\rm cheb}(y),
% \qquad
% y \in \mathbb{R}^m,
% where \Delta^{\rm cheb}: \mathbb{R}^m \to H is a linear map with an explicit expression given in the article. One considers a natural candidate for a functional estimating \gamma, namely the plug-in functional \delta^{\rm plug} := \gamma \circ \Delta^{\rm opti}: \mathbb{R}^m \to \mathbb{R}. Its error e(\delta^{\rm plug}) is computed via the second-order cone program stemming from  Proposition 4. Quite surprinsingly, it seems to agree with the minimal error e(\delta^{\rm opti}), meaning that, in this case, the plug-in functional is an optimal estimation functional, too.

C = U'*V;     % cross-gramian
D = U*(eye(m)-C*inv(C'*C)*C') + + V*inv(C'*C)*C';
b = W'*(g-D*U'*g);
Z = D'*W;
cvx_begin quiet
variable e_plug
variable sig(I_size,I_size) nonnegative
variable tau(I_size,I_size) nonnegative
minimize e_plug
subject to
sum(sig) == 1;
sum(tau) == 1;
for i=1:I_size
    V' * (W(:,i)-U*Z*sig(:,i)) == 0;
    V' * (U*Z(:,i)-W*tau(:,i)) == 0;
    epsilon*norm(W(:,i)-U*Z*sig(:,i)) <= ...
        e_plug + b'*sig(:,i) - g'*(W(:,i)-U*Z*sig(:,i));
    epsilon*norm(U*Z(:,i)-W*tau(:,i)) <= ...
        e_plug - b(i) - g'*(U*Z(:,i)-W*tau(:,i));
end
cvx_end
fprintf(['The error of the plug-in functional is %.3f, \n'...
    'which actually equals the minimal error of %.3f.'],...
    e_plug,e_opti)

%
%% References
% S. Foucart, "Optimal algorithms for nonlinear estimation with convex models", Preprint.
% M. Grant and S. Boyd, "CVX: Matlab software for disciplined convex programming", version 2.1. http://cvxr.com/cvx, March 2014.