<jupyter_start><jupyter_text>Computation of cutting planes: example 1  The set-up<jupyter_code>import numpy as np
import pandas as pd
import accpm
%load_ext autoreload
%autoreload 1
%aimport accpm<jupyter_output><empty_output><jupyter_text>$\DeclareMathOperator{\domain}{dom}\newcommand{\transpose}{\text{T}}\newcommand{\vec}[1]{\begin{pmatrix}1\end{pmatrix}}$  Example To test the computation of cutting planes we consider the unconstrained convex optimization problem\begin{align*}    &\text{minimize} \quad f_\text{obj}(x_0, x_1) = (x_0 - 5)^2 + (x_1 - 5)^2,\end{align*}and also the same problem with inequality constraints convex. That is, the problem\begin{align*}    &\text{minimize} \quad f_\text{obj}(x_0, x_1) = (x_0 - 5)^2 + (x_1 - 5)^2 \\    &\phantom{\text{minimize}} \quad f_0(x_0, x_1) =     a_0^\transpose x - b_0 = \vec{1\\0}^\transpose \vec{x_0\\x_1} - 20 = x_0 - 20 \leq 0\\    &\phantom{\text{minimize}} \quad f_1(x_0, x_1) =    a_1^\transpose x - b_1 = \vec{-1\\0}^\transpose \vec{x_0\\x_1} = -x_0 \leq 0\\    &\phantom{\text{minimize}} \quad f_2(x_0, x_1) =    a_2^\transpose x - b_2 = \vec{0\\1}^\transpose \vec{x_0\\x_1} - 20 = x_1 - 20 \leq 0 \\    &\phantom{\text{minimize}} \quad f_3(x_0, x_1) =     a_3^\transpose x - b_3 = \vec{0\\-1}^\transpose \vec{x_0\\x_1} = -x_1 \leq 0.\end{align*}In both cases it is clear that the solution is $x^\star = (x_1^\star, x_2^\star) = (5, 5)$. The ACCPM requires the gradients of the objective function and constraint functions, which are \begin{align*}    &\nabla f_\text{obj}(x_0, x_1) = \vec{2(x_0 - 5)\\2(x_1 - 5)}, \\    &\nabla f_0(x_0, x_1) = \vec{1\\0}, \quad \nabla f_1(x_0, x_1) = \vec{-1\\0}, \\    &\nabla f_2(x_0, x_1) = \vec{0\\1}, \quad \nabla f_3(x_0, x_1) = \vec{0\\-1}.\end{align*}We implement these functions as follows:<jupyter_code>def funcobj(x):
    return (x[0]-5)**2 + (x[1]-5)**2
    
def func0(x):
    return x[0] - 20

def func1(x):
    return -x[0]
    
def func2(x):
    return x[1] - 20
    
def func3(x):
    return -x[1]

def grad_funcobj(x):
    return np.array([2*(x[0] - 5), 2*(x[1] - 5)])

def grad_func0(x):
    return np.array([1, 0])

def grad_func1(x):
    return np.array([-1, 0])

def grad_func2(x):
    return np.array([0, 1])
    
def grad_func3(x):
    return np.array([0, -1])<jupyter_output><empty_output><jupyter_text>Here we analytically compute the initial few iterations for the unconstrained problem. The ACCPM requires that the initial polygon $\mathcal{P}_0$ (here I've abused terminology and by the initial polygon $\mathcal{P}_0$ I actually mean the system of linear inequalities $Ax \leq b$) contain at least some of the points we are interested in. For the purposes of this example we take\begin{align*}    A = \vec{a_0^\transpose\\a_1^\transpose\\a_2^\transpose\\a_3^\transpose}, b = \vec{20\\0\\20\\0}.\end{align*}Now, we start with $k=0$.  Now, $x^{(0)}_{ac}$ is the solution of the minimization problem \begin{equation*}    \min_{\domain \phi} \phi(x) = - \sum_{i=0}^{3}{\log{(b_i - a_i^\transpose x)}}.\end{equation*}So, we solve the problem\begin{align*}    &\phantom{iff}\nabla \phi(x) = \sum_{i=0}^{3    } \frac{1}{b_i - a_i^\transpose x}a_i = 0 \\    &\iff \frac{1}{20-x_0}\begin{bmatrix}1\\0\end{bmatrix} + \frac{1}{x_0}\begin{bmatrix}-1\\0\end{bmatrix} + \frac{1}{20-x_1}\begin{bmatrix}0\\1\end{bmatrix} + \frac{1}{x_1}\begin{bmatrix}0\\-1\end{bmatrix} = 0 \\    &\iff \frac{1}{20-x_0} - \frac{1}{x_0} = 0, \frac{1}{20-x_1} - \frac{1}{x_1} = 0 \\    &\iff x_0 = \frac{20}{2} = 10, x_1 = \frac{20}{2} = 10,\end{align*}and conclude $x^{(0)}_{ac} = (10, 10)$. We then query the oracle at $x^{(0)}_{ac}$. (Here, $f_\text{best} = f_\text{obj}(10, 10) = 50$ since this is the $0$-th iteration.) As there are no inequality constraints we have\begin{align*}    &a_4 = \nabla f_\text{obj}(10, 10) = \vec{10\\10}, \\    &b_4 = \nabla f_\text{obj}(10, 10)^\transpose \vec{10\\10} = \vec{10\\10}^\transpose \vec{10\\10} = 200,\end{align*}which we normalize to get\begin{align*}    &a_4 =  \frac{1}{\sqrt{100^2 + 100^2}} \nabla f_\text{obj}(10, 10)          = \vec{\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} } \approx \vec{0.7071 \\ 0.7071}, \\    &b_4 = \frac{1}{\sqrt{100^2 + 100^2}} \nabla f_\text{obj}(10, 10)^\transpose \vec{10\\10} = \vec{10\\10}^\transpose \vec{10\\10} = \frac{20}{\sqrt{2}} = 10\sqrt{2} \approx 14.1421,\end{align*}and therefore update\begin{align*}    A = \vec{a_0^\transpose\\a_1^\transpose\\a_2^\transpose\\a_3^\transpose\\ \frac{1}{\sqrt{2}} \;\; \frac{1}{\sqrt{2}}}, b = \vec{20\\0\\20\\0\\10\sqrt{2}}, k = 1.\end{align*}Now, $x^{(1)}_{ac}$ is the solution of the minimization problem \begin{equation*}    \min_{\domain \phi} \phi(x) = - \sum_{i=0}^{4}{\log{(b_i - a_i^\transpose x)}}.\end{equation*}So, we solve the problem\begin{align*}    &\phantom{iff}\nabla \phi(x) = \sum_{i=0}^{4    } \frac{1}{b_i - a_i^\transpose x}a_i = 0 \\    &\iff \frac{1}{20-x_0}\vec{1\\0} + \frac{1}{x_0}\vec{-1\\0} + \frac{1}{20-x_1}\vec{0\\1} + \frac{1}{x_1}\vec{0\\-1} + \frac{\sqrt{2}}{20-x_0-x_1} \vec{\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}} = 0\\    &\iff \frac{1}{20-x_0} - \frac{1}{x_0} + \frac{1}{20 - x_0- x_1}= 0, \frac{1}{20-x_1} - \frac{1}{x_1} + \frac{1}{20 - x_0- x_1} = 0 \\    &\iff x_0 = x_1 = 2(5 \pm \sqrt{5}) \approx 14.4721 \text{ or } 5.52786,\end{align*}and take $x^{(1)}_{ac} = (2(5-\sqrt{5}), 2(5-\sqrt{5})) \approx (5.52786, 5.52786)$. We then query theoracle at $x^{(1)}_{ac}$. Here $f_\text{obj}(x^{(1)}_{ac}) = f_\text{obj}(2(5-\sqrt{5}), 2(5-\sqrt{5})) = 90 - 40\sqrt{5} \approx 0.557281 \leq f_\text{best} = 50$ so we update $f_\text{best} = 90 - 40\sqrt{5} \approx 0.557281$ and therefore put (and normalize)  \begin{align*}    &a_5 = \frac{1}{\sqrt{2(10-4\sqrt{5})^2}} \nabla f_\text{obj}(x^{(1)}_{ac}) =            \frac{1}{\sqrt{2(10-4\sqrt{5})^2}} \nabla f_\text{obj}\vec{2(5-\sqrt{5}) \\ 2(5-\sqrt{5})} =           \frac{1}{\sqrt{2(10-4\sqrt{5})^2}} \vec{10-4\sqrt{5}\\10-4\sqrt{5}}           = \vec{\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}}\approx \vec{0.7071 \\ 0.7071}, \\    &b_5 = \frac{1}{\sqrt{2(10-4\sqrt{5})^2}} \nabla f_\text{obj}(x^{(1)}_{ac})^\transpose \vec{2(5-\sqrt{5}) \\ 2(5-\sqrt{5})} =            \frac{1}{\sqrt{2}} \vec{1\\1}^\transpose \vec{2(5-\sqrt{5}) \\ 2(5-\sqrt{5})} = 2\sqrt{2} (5 - \sqrt{5})           \approx 7.8176,\end{align*}updating $A$ and $b$ also. Iteration $k = 1$ is concluded by incrementing $k$. We see that this computation will only get more complicated. Therefore, in our implementation we test the computed analytic center by simply computing the norm of the gradient of the log barrier at the analytic center and see if it sufficiently small. We first test the unconstrained version of the problem.<jupyter_code>A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([20, 0, 20, 0])
accpm.accpm(A, b, funcobj, grad_funcobj, alpha=0.01, beta=0.7, 
            start=1, tol=10e-3, maxiter = 200, testing=1)<jupyter_output>-------- Starting ACCPM --------
Initially: b = [ 20.   0.  20.   0.] and A =
 [[ 1.  0.]
 [-1.  0.]
 [ 0.  1.]
 [ 0. -1.]]
--------------------------------
Entering iteration 0
At iteration 0 AC computation SUCCEEDED with AC [10. 10.] where
        a_cp = [ 0.71  0.71] and b_cp = [ 14.14]
--------------------------------
Entering iteration 1
At iteration 1 AC computation SUCCEEDED with AC [ 5.53  5.53] where
        a_cp = [ 0.71  0.71] and b_cp = [ 7.82]
--------------------------------
Entering iteration 2
At iteration 2 AC computation SUCCEEDED with AC [ 3.03  3.03] where
        a_cp = [-0.71 -0.71] and b_cp = [-5.58]
--------------------------------
Entering iteration 3
At iteration 3 AC computation FAILED with AC [ 2.26  5.72] where
        a_cp = [-0.97  0.25] and b_cp = [-2.06]
--------------------------------
Entering iteration 4
At iteration 4 AC computation FAILED with AC [ 2.43  7.31] where
        a_cp = [-0.74  0.67] and b_cp = [ 1.44]
--------------------------------
En[...]<jupyter_text>Next we test a trivial version of the inequality constrained problem where the inequality constraints and linear inequality constraints are the same<jupyter_code>A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([20, 0, 20, 0])
accpm.accpm(A, b, funcobj, grad_funcobj,
           (func0, func1, func2, func3), (grad_func0, grad_func1, grad_func2, grad_func3),
            alpha=0.01, beta=0.7, start=1, tol=10e-3, maxiter=200, testing=True)<jupyter_output>-------- Starting ACCPM --------
Initially: b = [ 20.   0.  20.   0.] and A =
 [[ 1.  0.]
 [-1.  0.]
 [ 0.  1.]
 [ 0. -1.]]
--------------------------------
Entering iteration 0
At iteration 0 AC computation SUCCEEDED with AC [10. 10.] where
        a_cp = [ 0.71  0.71] and b_cp = [ 14.14]
--------------------------------
Entering iteration 1
At iteration 1 AC computation SUCCEEDED with AC [ 5.53  5.53] where
        a_cp = [ 0.71  0.71] and b_cp = [ 7.82]
--------------------------------
Entering iteration 2
At iteration 2 AC computation SUCCEEDED with AC [ 3.03  3.03] where
        a_cp = [-0.71 -0.71] and b_cp = [-5.58]
--------------------------------
Entering iteration 3
At iteration 3 AC computation FAILED with AC [ 2.26  5.72] where
        a_cp = [-0.97  0.25] and b_cp = [-2.06]
--------------------------------
Entering iteration 4
At iteration 4 AC computation FAILED with AC [ 2.43  7.31] where
        a_cp = [-0.74  0.67] and b_cp = [ 1.44]
--------------------------------
En[...]<jupyter_text>Next we test a version of the inequality constrained problem where the initial polygon lies within the feasible region given by the inequality constraints.<jupyter_code>A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([5, 0, 5, 0])
accpm.accpm(A, b, funcobj, grad_funcobj,
           (func0, func1, func2, func3), (grad_func0, grad_func1, grad_func2, grad_func3),
            alpha=0.01, beta=0.7, start=1, tol=10e-3, maxiter=200, testing=True)<jupyter_output>-------- Starting ACCPM --------
Initially: b = [ 5.  0.  5.  0.] and A =
 [[ 1.  0.]
 [-1.  0.]
 [ 0.  1.]
 [ 0. -1.]]
--------------------------------
Entering iteration 0
At iteration 0 AC computation SUCCEEDED with AC [ 2.5  2.5] where
        a_cp = [-0.71 -0.71] and b_cp = [-3.54]
--------------------------------
Entering iteration 1
At iteration 1 AC computation SUCCEEDED with AC [ 3.62  3.62] where
        a_cp = [-0.71 -0.71] and b_cp = [-5.12]
--------------------------------
Entering iteration 2
At iteration 2 AC computation SUCCEEDED with AC [ 4.24  4.24] where
        a_cp = [-0.71 -0.71] and b_cp = [-6.]
--------------------------------
Entering iteration 3
At iteration 3 AC computation SUCCEEDED with AC [ 4.59  4.59] where
        a_cp = [-0.71 -0.71] and b_cp = [-6.49]
--------------------------------
Entering iteration 4
At iteration 4 AC computation SUCCEEDED with AC [ 4.78  4.78] where
        a_cp = [-0.71 -0.71] and b_cp = [-6.75]
--------------------------------
E[...]<jupyter_text>We now test a version of the inequality constrained problem where the initial polyhedron contains the feasible region but also regions that are infeasible.<jupyter_code>A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([30, 0, 30, 0])
accpm.accpm(A, b, funcobj, grad_funcobj,
           (func0, func1, func2, func3), (grad_func0, grad_func1, grad_func2, grad_func3),
            alpha=0.01, beta=0.7, start=1, tol=10e-3, maxiter=200, testing=True)<jupyter_output>-------- Starting ACCPM --------
Initially: b = [ 30.   0.  30.   0.] and A =
 [[ 1.  0.]
 [-1.  0.]
 [ 0.  1.]
 [ 0. -1.]]
--------------------------------
Entering iteration 0
At iteration 0 AC computation SUCCEEDED with AC [ 15.  15.] where
        a_cp = [ 0.71  0.71] and b_cp = [ 21.21]
--------------------------------
Entering iteration 1
At iteration 1 AC computation SUCCEEDED with AC [ 8.29  8.29] where
        a_cp = [ 0.71  0.71] and b_cp = [ 11.72]
--------------------------------
Entering iteration 2
At iteration 2 AC computation SUCCEEDED with AC [ 4.54  4.54] where
        a_cp = [-0.71 -0.71] and b_cp = [-6.42]
--------------------------------
Entering iteration 3
At iteration 3 AC computation FAILED with AC [ 3.21  8.02] where
        a_cp = [-0.51  0.86] and b_cp = [ 3.58]
--------------------------------
Entering iteration 4
At iteration 4 AC computation SUCCEEDED with AC [ 10.15   3.18] where
        a_cp = [ 0.94 -0.33] and b_cp = [ 5.82]
---------------------------[...]<jupyter_text>We test a version of the inequality constrained problem where the initial polyhedron is moderately large.<jupyter_code>A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([100, 0, 100, 0])
accpm.accpm(A, b, funcobj, grad_funcobj,
           (func0, func1, func2, func3), (grad_func0, grad_func1, grad_func2, grad_func3),
            alpha=0.01, beta=0.7, start=1, tol=10e-3, maxiter=200, testing=True)<jupyter_output>-------- Starting ACCPM --------
Initially: b = [ 100.    0.  100.    0.] and A =
 [[ 1.  0.]
 [-1.  0.]
 [ 0.  1.]
 [ 0. -1.]]
--------------------------------
Entering iteration 0
At iteration 0 AC computation SUCCEEDED with AC [ 50.  50.] where
        a_cp = [ 1.  0.] and b_cp = [ 20.]
--------------------------------
Entering iteration 1
At iteration 1 AC computation SUCCEEDED with AC [  9.45  50.  ] where
        a_cp = [ 0.  1.] and b_cp = [ 20.]
--------------------------------
Entering iteration 2
At iteration 2 AC computation SUCCEEDED with AC [ 9.45  9.45] where
        a_cp = [ 0.71  0.71] and b_cp = [ 13.36]
--------------------------------
Entering iteration 3
At iteration 3 AC computation FAILED with AC [ 13.48   6.35] where
        a_cp = [ 0.99  0.16] and b_cp = [ 12.32]
--------------------------------
Entering iteration 4
At iteration 4 AC computation FAILED with AC [ 13.51   6.27] where
        a_cp = [ 0.99  0.15] and b_cp = [ 12.29]
-------------------------------[...]<jupyter_text>Finally we test a version of the inequality constrained problem where the initial polyhedron is very large. We observe that the algorithm does not converge. It is possible more iterations, or pruning of the redundant inequalities, would allow for convergence, but in practice, it is unlikely this situation will arise.<jupyter_code>A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([1000, 0, 1000, 0])
accpm.accpm(A, b, funcobj, grad_funcobj,
           (func0, func1, func2, func3), (grad_func0, grad_func1, grad_func2, grad_func3),
            alpha=0.01, beta=0.7, start=1, tol=10e-3, maxiter=200, testing=True)
A = np.array([[1, 0],[-1,0],[0,1],[0,-1]])
b = np.array([1000, 0, 1000, 0])
accpm.accpm(A, b, funcobj, grad_funcobj,
           (func0, func1, func2, func3), (grad_func0, grad_func1, grad_func2, grad_func3),
            alpha=0.01, beta=0.7, start=1, tol=10e-3, maxiter=500, testing=True)<jupyter_output>-------- Starting ACCPM --------
Initially: b = [ 1000.     0.  1000.     0.] and A =
 [[ 1.  0.]
 [-1.  0.]
 [ 0.  1.]
 [ 0. -1.]]
--------------------------------
Entering iteration 0
At iteration 0 AC computation SUCCEEDED with AC [ 499.96  499.96] where
        a_cp = [ 1.  0.] and b_cp = [ 20.]
--------------------------------
Entering iteration 1
At iteration 1 AC computation SUCCEEDED with AC [   9.95  500.  ] where
        a_cp = [ 0.  1.] and b_cp = [ 20.]
--------------------------------
Entering iteration 2
At iteration 2 AC computation SUCCEEDED with AC [ 9.95  9.95] where
        a_cp = [ 0.71  0.71] and b_cp = [ 14.07]
--------------------------------
Entering iteration 3
At iteration 3 AC computation FAILED with AC [ 11.95   9.2 ] where
        a_cp = [ 0.86  0.52] and b_cp = [ 13.94]
--------------------------------
Entering iteration 4
At iteration 4 AC computation FAILED with AC [ 14.3    6.32] where
        a_cp = [ 0.99  0.14] and b_cp = [ 12.96]
-------------------[...]
{'fraction_non_alphanumeric': 0.2483361855866134, 'fraction_numerical': 0.08689864993344742, 'mean_word_length': 2.2578979971092297, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '509937', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6515, 'n_tokens_neox': 5822, 'n_words': 1620, 'path': 'projects/david/lab/test_accpm_toy.ipynb', 'repo_name': 'chengsoonong/mclass-sky', 'size': '15777'}