Clone wiki

EmuMore / shape_constraints_polynomial

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);
    # The integral could be a non-adaptive quadrature
    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
upstream network with inlets

Updated