Using Mathematica to test and implement a new algorithm for inequality and equality constrained optimization problems

Philibert Kongtcheu

PFK Technologies

Problem Statement:
(1){({{Find a Global Minimum f((Overscript[x, ∼]))at (Overscript[x, ∼]) for f:Ω⊂(R^n)→R subject to the conditions}, {g(x)≤(Subscript[0, m]),g:Ω⊂(R^n)→(R^m)}, {h(x)=(Subscript[0, r]),h:Ω⊂(R^n)→(R^r)}, {and f,g, h are twice continuously differentiable }})

Or

(1'){({{Find a Global Minimum f((Overscript[x, ∼]))at (Overscript[x, ∼]) for f:Ω⊂(R^n)→R subject to the conditions}, {(Subscript[L, m])≤g(x)≤(Subscript[H, m]),g:Ω⊂(R^n)→(R^m),{(Subscript[L, m]),(Subscript[H, m])}⊂(R^m)}, {h(x)=(Subscript[0, r]),h:Ω⊂(R^n)→(R^r)}, {and f,g, h are twice continuously differentiable }})

Which can also be written as:

(1''){({{Find a Global Minimum f((Overscript[x, ∼]))at (Overscript[x, ∼]) for f:Ω⊂(R^n)→R subject to the conditions}, {(Subscript[L, m])≤g(x)≤(Subscript[H, m]),g:Ω⊂(R^n)→(R^m),{(Subscript[L, m]),(Subscript[H, m])}⊂(R^m)}, {and f,g,are twice continuously differentiable }})

Toobox Tool 1: The Karush Kuhn Tucker (KKT) Theorem

Toobox Tool 1: The KKT - Necessary Conditions

(2){({{For a point (Overscript[x, ∼])to yield a minimum at f((Overscript[x, ∼]))at (Overscript[x, ∼]) for f:Ω⊂(R^n)→R subject to the conditions}, {g(x)≤(Subscript[0, m]),g:Ω⊂(R^n)→(R^m)}, {h(x)=(Subscript[0, r]),h:Ω⊂(R^n)→(R^r)}, {and f,g, h are twice continuously differentiable,the following conditions must hold: }})

Alternative KKT Formulation:

(2'){({{For a point (Overscript[x, ∼])to yield a minimum at f((Overscript[x, ∼]))at (Overscript[x, ∼]) for f:Ω⊂(R^n)→R subject to the conditions}, {(Subscript[L, m])≤g(x)≤(Subscript[H, m]),g:Ω⊂(R^n)→(R^m),{(Subscript[L, m]),(Subscript[H, m])}⊂(R^m)}, {h(x)=(Subscript[0, r]),h:Ω⊂(R^n)→(R^r)}, {and f,g, h are twice continuously differentiable,the following conditions must hold: }})

1. Stationarity;As a linear application,there exists λ in (R^r), such that:
dKKT=D(f+<λ,h>+<(SuperPlus[μ]),g-(Subscript[H, m])>-<(μ^−),g-(Subscript[L, m])>)((Overscript[x, ∼]))=0;
We may therefore subsequently note μ=(SuperPlus[μ])-(SuperMinus[μ]);
2. Primal feasibility Equality constraints conditions:h((Overscript[x, ∼]))=0;
3. Primal feasibility inequality constraints conditions:|g((Overscript[x, ˜]))-(1/2)((Subscript[H, m])+(Subscript[L, m]))|≤(1/2)((Subscript[H, m])-(Subscript[L, m]));
4. Complementary slackness or switching conditions:(SuperPlus[μ])*(g((Overscript[x, ˜]))-(Subscript[H, m]))=0;(μ^−)*(g((Overscript[x, ˜]))-(Subscript[L, m]))=0;
5. Dual Feasibility sign of inequality multipliers :(SuperPlus[μ]),(SuperMinus[μ])in  (Subscript[R, +]^m) .

Toobox Tool 2: The Second Order Necessary or Sufficient Conditions

Theorem: Second Order Necessary Conditions

We note : (Subscript[N, ϵ]) (x)   (Overscript[=, Δ])NullSpace[(({{(Subscript[(((∂(Subscript[h, i]))/(∂(Subscript[x, j])))), 1≤i≤r,1≤j≤n])}, {(Subscript[(((∂(Subscript[g, i']))/(∂(Subscript[x, j])))), i' in active contraints set ϵ  ,1≤j≤n])}}))];
For a point (SuperStar[x])*such that ((SuperStar[x]),(SuperStar[μ]),(SuperStar[λ])) satisfy the necessary KKT conditions to be a local minimum of f and where f,g,h are continuously twice differentiable at (SuperStar[x]),the following condition is a necessary condition for (SuperStar[x])to be a local minimum of f subject to the constraints (Subscript[h, i])((SuperStar[x]))=0,i=1,2,…,p and (Subscript[g, j])((SuperStar[x]))≤0,j=1,2,…,m:
Rank[(Subscript[N, ϵ]) ((SuperStar[x])) ]=0  or  as a symmetric matrix,(Subscript[N, ϵ]) ((SuperStar[x]))(D^2)(f+<(SuperStar[μ]),g>+<(SuperStar[λ]),h>)(Subscript[N, ϵ]) ((((SuperStar[x])))^T)is positive semi-definite.


This theorem helps give us tools to broaden our investigations beyond convex optimization problems and not limit ourselves to so-called invex functions.

Theorem:Second Order Sufficient Conditions

For a point (SuperStar[x])*such that ((SuperStar[x]),(SuperStar[μ]),(SuperStar[λ])) satisfy the necessary KKT conditions to be a local minimum of f and where f,g,h are continuously twice differentiable at (SuperStar[x]),the following condition is a sufficient condition for (SuperStar[x])to be a local minimum of f subject to the constraints (Subscript[h, i])((SuperStar[x]))=0,i=1,2,…,p and (Subscript[g, j])((SuperStar[x]))≤0,j=1,2,…,m:
1.Rank[(Subscript[N, ϵ]) ((SuperStar[x])) ]>0;  
2. as a symmetric matrix,(Subscript[N, ϵ]) ((SuperStar[x]))(D^2)(f+<(SuperStar[μ]),g>+<(SuperStar[λ]),h>)(Subscript[N, ϵ]) ((((SuperStar[x])))^T)is positive definite and;
3.For any j such that (Subscript[g, j])((SuperStar[x]))=0,(Subscript[μ, j])>0.

Mathematica Issue: No Seamless Computation of an analytic NullSpace function

Toolbox Issues: Computational Tractability - Curse of Dimensionality - N vs NP

•	We note s the number of coordinates j such that (Subscript[μ, j])>0. s varies from 0 to m. 
•	We note (Subscript[ε, s]) an m-dimensional vector made of 0s and 1s such that the sum of its non zero terms equals s. The position j of each number 1 in this system indicates a j such that (Subscript[g, j])(x*)=0. For other coordinates, λj=0. There are such terms (Subscript[ε, s]). 
The resolution of KKT equations and inequalities as formulated hides intractable complexity. As we shall see, computational complexity grows exponentially with m. The first order KKT conditions potentially outline a system of:
•	n+r+m  equations and n+p+m variables;
•	2m inequalities
We can break these conditions along the number of indices for which inequalities reduce to equalities. s (or alternative called pp) is the number of indices j where the number of coordinates j such that μj>0 and leads to the inequality (Subscript[g, j])≤0 reducing to (Subscript[g, j])=0. The KKT conditions are transformed into 2m systems of equations and inequalities sub-divided for any fixed s ranging from zero to m into  (({{m}, {s}}))=((m!)/(s!(m−s)!))systems of equations and inequalities each representing:
•	n+p+s equations  and n+p+s variables;
•	m inequalities; 
Each system of  n+r+s equations with n+r+s variables is then solved and one checks that the solutions verify the m inequalities.

A.  Serial implementation complexity

The resolution of each such system using the best equation solvers generally proceeds in O((((n+p+s))^3)) time. Therefore the total computational time is (∑_(s=0)^m)(({{m}, {s}}))((O(n+r+s))^3).
                                             But  (∑_(s=0)^m)(({{m}, {s}}))((O(n+r+s))^3)≈(2^m)O((((n+r))^3))

N vs NP problem...Hopeless...

Prior Art: Interior Point Methods

Overview

Interior point methods (also referred to as barrier methods) are a certain class of algorithms to solve linear and nonlinear convex optimization problems.

The interior point method was invented by John von Neumann.Von Neumann suggested a new method of linear programming,using the homogeneous linear system of Gordan (1873) which was later popularized by Karmarkar's algorithm,developed by Narendra Karmarkar in 1984 for linear programming.The method consists of a self-concordant barrier function used to encode the convex set.Contrary to the simplex method,it reaches an optimal solution by traversing the interior of the feasible region.Any convex optimization problem can be transformed into minimizing (or maximizing) a linear function over a convex set.The idea of encoding the feasible set using a barrier and designing barrier methods was studied in the early 1960s by,amongst others,Anthony V.Fiacco and Garth P.McCormick.These ideas were mainly developed for general nonlinear programming,but they were later abandoned due to the presence of more competitive methods for this class of problems (e.g.sequential quadratic programming).

Our Approach

Example 1: Quadratic Minimization

Problem Statement

This problem can be generally stated analytically as:Find the solution to the constrained quadratic optimization problem
({{Min (1/2)(W^t)ΣW +(C^t)Wsubject  to constraints }, {{({{AW−B=0;}, {|W−(((Overscript[M, ¯])+(Underscript[M, ¯]))/2)|≤(((Overscript[M, ¯])−(Underscript[M, ¯]))/2)}})}})
with  (Underscript[M, ¯]),W,C,(Overscript[M, ¯])∈(R^m),(Underscript[M, ¯])≤W≤(Overscript[M, ¯]),A∈(R^s)×(R^m);Rank(A)=s≤m,B∈(R^s)

Analytic Solution Derivation

Local Extrema Derivation

(Subscript[W, q])=−(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));     (Subscript[γ, p])=−(Subscript[s, pp])(Subscript[M, p])−(Subscript[A, p]^t)λ'−(Subscript[C, p])+(Subscript[s, pq])(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));

Where p+q=m; 
(Subscript[W, q])is the vector of non active constraints and 
(Subscript[γ, p]) represents the slacking parameters vector for the active constraints. 

This provides the general expression of the KKT local extremum.In order for them to be admissible,they have to verify the admissibility conditions.

Solution

Local Extrema Derivation

(Subscript[W, q])=−(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));     (Subscript[γ, p])=−(Subscript[s, pp])(Subscript[M, p])−(Subscript[A, p]^t)λ'−(Subscript[C, p])+(Subscript[s, pq])(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));

Where p+q=m; 
Wq is the vector of non active constraints and 
(Subscript[γ, p]) represents the slacking parameters vector for the active constraints. 

This provides the general expression of the KKT local minima.In order for them to be admissible,they have to verify the admissibility conditions.

Case Study 1

Solution

Local Extrema Derivation

(Subscript[W, q])=−(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));     (Subscript[γ, p])=−(Subscript[s, pp])(Subscript[M, p])−(Subscript[A, p]^t)λ'−(Subscript[C, p])+(Subscript[s, pq])(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));

Where p+q=m; 
Wq is the vector of non active constraints and 
(Subscript[γ, p]) represents the slacking parameters vector for the active constraints. 

This provides the general expression of the KKT local minima.In order for them to be admissible,they have to verify the admissibility conditions.

Case Study

This example is a classic portfolio optimization problem where one seeks to minimize a portfolio’s variance subject to allocations size constraints and a target return.The data comes from Stock quotes of a 5 large Caps’ 3-month Data July 30 to Sept 30,2007 .

We use the covariance matrix Sig extracted from above,and the individual return vector R.The components in the allocations vector W must range between-1 and 1 and sum to 1 (-1≤Wi≤1,Σ Wi=1),the target return is 0.674% (0.006739289957201789`) and we wish to minimize the variance WtSigW.
5=m is the number of assets.

Results & Extensions: See  Corresponding Notebook or (kongtcheu@pfktech.com)

Example 2: Linear Maximization with Quadractic Constraints

Problem Statement

This problem can be generally stated analytically as:Find the solution to the constrained quadratic optimization problem
({{Max (R^t)Wsubject  to constraints }, {{({{AW−B=0;(W^t)ΣW =(σ^2);}, {|W−(((Overscript[M, ¯])+(Underscript[M, ¯]))/2)|≤(((Overscript[M, ¯])−(Underscript[M, ¯]))/2)}})}})
with  (Underscript[M, ¯]),W,C,(Overscript[M, ¯])∈(R^m),(Underscript[M, ¯])≤W≤(Overscript[M, ¯]),A∈(R^s)×(R^m);Rank(A)=s≤m,B∈(R^s)

Analytic Solution Derivation

Local Extrema Derivation

(Subscript[W, q])=−(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));     (Subscript[γ, p])=−(Subscript[s, pp])(Subscript[M, p])−(Subscript[A, p]^t)λ'−(Subscript[C, p])+(Subscript[s, pq])(Subscript[s, qq]^(-1))((Subscript[s, pq]^t)(Subscript[M, p])+(Subscript[A, q]^t)λ'+(Subscript[C, q]));

Where p+q=m; 
(Subscript[W, q])is the vector of non active constraints and 
(Subscript[γ, p]) represents the slacking parameters vector for the active constraints. 

This provides the general expression of the KKT local extremum.In order for them to be admissible,they have to verify the admissibility conditions.

Solution

Local Extrema Derivation

W=((((a^t)Σb)/((a^t)Σa))−(Sqrt[(1/((a^t)Σa))((σ^2)−((((a^t)Σa)((b^t)Σb)−((((a^t)Σb ))^2))/((a^t)Σa)))])) a−b;({{(Subscript[γ, p])=λ(((Subscript[s, pq])(Subscript[s, qq]^(−1))(Subscript[s, pq]^t)−(Subscript[s, pp])+((Subscript[s, pq])(Subscript[s, qq]^(−1))(Subscript[A, q]^t)−(Subscript[A, p]^t))(Subscript[χ, s]^(−1))((Subscript[A, p])−(Subscript[A, q])(Subscript[s, qq]^(−1))(Subscript[s, pq]^t)))(Subscript[M, p])−((Subscript[s, pq])(Subscript[s, qq]^(−1))(Subscript[A, q]^t)−(Subscript[A, p]^t))(Subscript[χ, s]^(−1))B)}, {+((Subscript[s, pq])(Subscript[s, qq]^(−1))(Subscript[R, q])−(Subscript[R, p]))−((Subscript[s, pq])(Subscript[s, qq]^(−1))(Subscript[A, q]^t)−(Subscript[A, p]^t))(Subscript[χ, s]^(−1))(Subscript[A, q])(Subscript[s, qq]^(−1))(Subscript[R, q])}})




Where p+q=m; 
Wq is the vector of non active constraints and 
(Subscript[γ, p]) represents the slacking parameters vector for the active constraints. 
(Subscript[χ, s])=(Subscript[A, q])(Subscript[s, qq]^(−1))(Subscript[A, q]^t);a=(({{0}, {−(Subscript[s, qq]^(−1))((Subscript[I, q])−(Subscript[A, q]^t)(Subscript[χ, s]^(−1))(Subscript[A, q])(Subscript[s, qq]^(−1)))(Subscript[R, q])}}));b=(({{−(Subscript[I, p])}, {(Subscript[s, qq]^(−1))((Subscript[s, pq]^t)+(Subscript[A, q]^t)(Subscript[χ, s]^(−1))((Subscript[A, p])−(Subscript[A, q])(Subscript[s, qq]^(−1))(Subscript[s, pq]^t)))}}))(Subscript[M, p])−(({{(Subscript[0, p])}, {(Subscript[s, qq]^(−1))(Subscript[A, q]^t)(Subscript[χ, s]^(−1))B}}));
(1/λ)=(((a^t)Σb)/((a^t)Σa))-(Sqrt[(1/((a^t)Σa))((σ^2)−((((a^t)Σa)((b^t)Σb)−((((a^t)Σb ))^2))/((a^t)Σa)))])

This provides the general expression of the KKT local minima.In order for them to be admissible,they have to verify the admissibility conditions.

Investigating Convergence Issues

Open Question

What is the natural boundary of the number of cases that need be be sample to obtain a global optimum?
A hard Math Estimation may not be easy
Implementation probing +systematic random case generation

Problem Statement

Before KKT equations and inequalities are solved:

Dimension Reduction from permutations invariance or other dimension reducing transformations

Change of variable to facilitate treatment of inequality constraints or to facilitate deducing if the Hessian is positive definite (for a minimum) or negative definite (for a maximum).

In addition:Use the Order Relationship <<to Eliminate unnecessary cases from evaluation list

Matrix Inverse Computational Reduction By exploitation of Relationships on the Order Relationship <<dependency tree

Grouping of Cases i.e.along same inverse matrix requirements or other parameters