# Shape constraints for global polynomials

Given a domain $$\Omega \in \mathbb{R}^d$$, define a finite partition with $$K$$ subdomains

\begin{equation*} P_K(\Omega) = \left\lbrace\Omega_k : \Omega_k \cap\Omega_{k'} = \emptyset, k\neq k' \wedge {\scriptsize\bigcup}_k \Omega_k = \Omega\right\rbrace \end{equation*}

Find $$p(x) \in \Pi(\Omega)$$ that minimizes

\begin{equation*} \operatorname{Loss}\left(\lbrace y_i\rbrace, p[\lbrace x_i\rbrace]\right) \end{equation*}

subject to the functional inequality constraints

\begin{equation*} h_k^s(\Omega_k, p) := \int_{\Omega_k}\left( p^{(s)}(x) - \vert p^{(s)}(x) \vert \right) dx \geq 0 \end{equation*}

## A simple 1D example

Definition of integrand in functional $$h_k^s$$

function v = h_integrand (s, p, t)
# derivate the polynomial |s| times
for i = 1:abs(s)
p = polyder (p);
endfor
# evaluate the polynomial
p = polyval (p, t) * sign (s);
# integrand
v = p - abs (p);
endfunction


Definition of functional

function v = h (s, Omega, p)
func = @(t) h_integrand (s, p, t);
v    = quadgk (func, Omega(1), Omega(2));
endfunction


Generate data for the example

n = 100;
t = linspace (-1, 1, n).';
y = tanh (2 * t) + 0.1 * randn (n, 1);


Define the optimization problem and solve

# The parameter vector is composed of D+1 elements describing a polynomial of
# degree D followed by 1 element describing the position of the inflection point.
#
D  = 10; % polynomial degree
O  = D + 1;
p0 = polyfit (t, y, D);
p0(end+1) = 2 * rand - 1; % guess of inflection point
p0 = p0.';

# Data for the gradient of cost
Vd = vander (t, O);
K = Vd.' * Vd;
Vy = Vd.' * y;
H = blkdiag (K, 1);
cost ={@(x) x(1:O).'* ( K * x(1:O) - 2 * Vy );
@(x) [2*(K*x(1:O)-Vy); 0];
@(x) H};
# inequalities
ineq =@(x) [ h(2, [-1 x(end)], x(1:O)); ... % concave 1st half
h(-2, [x(end) 1], x(1:O))];    % convex  2nd half

# Solve using SQP, I believe that with a fixed quadrature this becomes QP
tic
[p, obj, infocode] = sqp (p0, cost, [], ineq);
toc
t_inf = p(end);   % inflection point
p     = p(1:O); % polynomial coefficients

# Plot solution
hp = plot (t, y, 'x', ...
t, polyval (p, t), '-', ...
t_inf, polyval (p, t_inf), 'og');
set (hp(end), 'markerfacecolor', 'auto');
axis tight


Updated