<jupyter_start><jupyter_text>Constrained optimization: equality constraints By <NAME> and Johnny  Sequential Quadratic Programming<jupyter_code>import numpy as np
def f(x):
    return np.array([np.power(np.e, 3 * x[0]) + np.power(np.e, -4 * x[1])], dtype=np.float64)


def h(x):
    return np.array([np.power(x[0], 2) + np.power(x[1], 2) - 1], dtype=np.float64)
def df(x):
    d1 = 3 * np.power(np.e, 3 * x[0])
    d2 = -4 * np.power(np.e, -4 * x[1])
    return np.array([d1, d2])


def dh(x):
    d1 = 2 * x[0]
    d2 = 2 * x[1]
    return np.array([d1, d2])


def dl(x, lamb):
    return df(x) - lamb * dh(x)
def d2f(x):
    d11 = 9 * np.power(np.e, 3 * x[0])
    d12 = 0
    d21 = 0
    d22 = 16 * np.power(np.e, -4 * x[1])
    return np.array([[d11, d12], [d21, d22]])


def d2h(x):
    d11 = 2
    d12 = 0
    d21 = 0
    d22 = 2
    return np.array([[d11, d12], [d21, d22]])


def d2l(x):
    return d2f(x) - lamb * d2h(x)
def d2la(x, lamb):
    G = d2f(x) - lamb * d2h(x)
    C = -dh(x)

    matrix = np.zeros((3, 3))
    matrix[:2, :2] = G
    matrix[-1, :2] = C
    matrix[:2, -1] = C
    
    return matrix<jupyter_output><empty_output><jupyter_text>1. The stopping criterion we have estalished is to check that the norm of $\nabla xL$ is major than epsilon.<jupyter_code>def solve(x, lamb, max_iter=1000, eps = 10e-16):
    i = 0

    dlx = dl(x, lamb)

    while i < max_iter and np.linalg.norm(dlx) > eps:
        A = d2la(x, lamb)
        dlx = dl(x, lamb)
        b = np.concatenate((-dlx, h(x)))
        d = np.linalg.solve(A, b)

        x += d[:-1]
        lamb += d[-1]
        i += 1

    return x, lamb
x = np.array([-1., 1.])

lamb = -1

print(solve(x, lamb))<jupyter_output>(array([-0.74833549,  0.66332043]), -0.21232493554997134)
(array([-0.74833549,  0.66332043]), -0.2123249355499713)<jupyter_text>1. We can can see that the algorithm with initial point closed to the minimum can be converged.<jupyter_code>def merit(x, p=10):
    return f(x) + p * np.power(h(x), 2)

def dmerit(x, p=10):
    d1 = df(x)
    d2 = 2 * p * dh(x)
    return d1 + d2<jupyter_output><empty_output><jupyter_text>2. Setting x = (100, 100).<jupyter_code>x = np.array([100., 100.])

lamb = -1

print(solve(x, lamb))<jupyter_output>(array([nan, nan]), nan)<jupyter_text>2. We can see that with initial point far away from minimum pont cannot be solved by Newton-like iterarion.  3. Implement gradient descent method with Merit function.<jupyter_code>def gradient_descend_2d(dmerit, x, alpha=1, eps=1e-6, max_iter=1000000):
    i = 0

    while i < max_iter:
        dm = dmerit(x)
        dm_norm = np.linalg.norm(dm)
        dm_normalized = dm / dm_norm

        x -= alpha * dm_normalized
        i += 1

    return x<jupyter_output><empty_output><jupyter_text>In order to find the minimum, we have used the gradient descent method. In this gradient descent methods, we have penalized the constraint function quadratically when the initial point is far away from the minimum.  4. Gradient descent and Newton-like iteration with Merit function.<jupyter_code>x = np.array([100., 100.])

minimum = gradient_descend_2d(dmerit, x, max_iter=100)

print(minimum)

lamb = -1

print(solve(minimum, lamb))<jupyter_output>[ 1.55672317 97.89137348]
(array([-0.74833549,  0.66332043]), -0.2123249355499713)
{'fraction_non_alphanumeric': 0.15629742033383914, 'fraction_numerical': 0.08679817905918058, 'mean_word_length': 3.054120541205412, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 24, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '411135', 'license': 'MIT', 'n_tokens_mistral': 1550, 'n_tokens_neox': 1369, 'n_words': 357, 'path': 'OPT/Lab/Lab4_QiJun_Johnny.ipynb', 'repo_name': 'johnnync13/DataScienceUB', 'size': '3298'}