# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # CVXOPT

# ## Task 1
# Machine Learning tasks are typically thought of optimization problems, e.g. minimizing an error function or maximizing a probability. Ideally, the optimization problem turns out to be convex, which implies that any local minimum is the global minimum of the formulation, and what is even more important, we can.  In the following, it will be assumed that you have some basic knowledge about convex optimization. The intention of this task is to familiarize ourselves with CVXOPT, one of the most-widely used convex optimization toolboxes.

# a)  Go to `cvxopt.org` and follow the installation instructions for your distribution. For conda, you need to run
# `conda install -c conda-forge cvxopt`

# b) Skim through the **Examples** section on `cvxopt.org` to get an overview of the functionality of the different solvers of CVXOPT.

from cvxopt import matrix, solvers
import numpy as np


# c) Implement a function `minsq` which expects a NumPy array `A` of shape `(m,n)` and a NumPy array `y` of shape `(m,)` as its arguments and returns a NumPy array `x` of shape `(n,)` that solves the following problem.
#
# <center>$\mathrm{min_\mathbf{x}} \|\mathbf A\mathbf{x}-\mathbf{y}\|$.</center>
#
# Test your function by feeding it with appropriate inputs and comparing the results with the ones you get by using `np.linalg.pinv`. Experiment by adding white Gaussian noise to `y`. If CVXOPT does not accept your NumPy arrays, try casting them to `double`.

# +
def minsq(A, y):
    P=matrix(np.dot(A.T,A).astype('double'))
    q=matrix(-np.dot(A.T,y).astype('double'))
    x=solvers.qp(P,q)
    return np.array(x['x'])

A=np.array([[10, 40],[20, 0],[-30, 40]])
y=np.array([50,20,10])+np.random.randn(3,)

print('A:', A)
print('y:', y)
print('x:', minsq(A,y).squeeze())
print('np.dot(pinv(A),y):', np.dot(np.linalg.pinv(A),y))


# -

# d) Consider the equation (8.30) in the lecture notes. Implement a function `solvedualsvm(H,y)` that returns the solution `lambda_star` of the dual SVM problem by means of CVXOPT.

def solvedualsvm(H,y):
    y=y.squeeze()
    n=len(y)
    G=-np.eye(n).astype('double')
    A=y.reshape(1,n).astype('double')
    h=np.zeros((n,)).astype('double')
    b=np.zeros((1,)).astype('double')
    P=H.astype('double')
    q=-np.ones((n,)).astype('double')
    lambda_star=solvers.qp(matrix(P),matrix(q),matrix(G),matrix(h),matrix(A),matrix(b))
    return lambda_star['x']


# Test your function with the training data
# 		\begin{equation*}
# 		\begin{split}
# 		\mathbf{x}_1=\begin{bmatrix}-1\\-1\end{bmatrix},y_1=-1,&\ \mathbf{x}_2=\begin{bmatrix}-2\\-2\end{bmatrix},y_2=-1,\\
# 		\mathbf{x}_3=\begin{bmatrix}1\\1\end{bmatrix},y_3=1,&\ \mathbf{x}_4=\begin{bmatrix}2\\2\end{bmatrix},y_4=1,\
# 		\end{split}.
# 		\end{equation*}
# 		Verify that the KKT conditions with respect to the support vectors are in line with what you expect. In the next lab course, we will use this function to implement linear and kernel SVM.

X=np.array([[-1,-2,1,2],[-1,-2,1,2]])
y=np.array([-1,-1,1,1])
H=np.dot(np.dot(np.dot(np.diag(y),X.T),X),np.diag(y))
lambda_star=solvedualsvm(H,y)
print(lambda_star)

# Only the KKT coefficients that belong to the support vectors are significantly larger than 0.

StarcoderdataJupyterScriptsDedupFiltered