pyrobolearn/tutorials/math/PCA.ipynb

PCA

In this notebook, I will review a little bit about PCA, implement recursive PCA, how PCA can be viewed as an optimization problem, and implement a constrained version of this optimization process for PCA applied on time-series by including a roughness penalty.

Generalized PCA will not be considered after careful consideration. Note: "Generalized PCA" is "an algebro-geometric solution to the problem of segmenting an unknown number of subspaces of unknown and varying dimensions from sample data points".

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d

from sklearn.decomposition import PCA
%matplotlib inline

In [ ]:

class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)

Let's first start to apply PCA on the following toy example.

In [ ]:

# Toy data
X = np.array([[-1, -1],
              [-2, -1],
              [-3, -2],
              [1, 1],
              [2, 1],
              [3, 2]], dtype=np.float64) # NxT

# Plot
plt.figure(figsize=(5,4))
plt.scatter(X[:,0], X[:,1], color='b')
plt.arrow(0, 0, 1, 0, length_includes_head = True, head_width = 0.15, color='k')
plt.arrow(0, 0, 0, 1, length_includes_head = True, head_width = 0.15, color='')
plt.show()

In [ ]:

### PCA
n_components = 2 # number of principal axes that we want to keep

## Let's compute PCA from scratch
# 1. Center the data
mean = X.mean(axis=0)
X -= mean # note that the data was already centered
N = X.shape[0]

# 2. Compute the covariance matrix
CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))

# 3. Compute the eigenvectors of this covariance matrix
# np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix
evals, evecs = np.linalg.eigh(CovX)

# 4. Sort the eigenvalues (in decreasing order) and eigenvectors
idx = np.argsort(evals)[::-1]
evals = evals[idx]
evecs = evecs[:,idx]

# 5. Form the projection matrix
P = evecs[:,:n_components]
print(P)

# 5. Project the data using the projection matrix
# This is the same as rotating the matrix X using P
Y = X.dot(P)
    
# 6. Compare it with standard PCA
pca = PCA(n_components=n_components)
pca = pca.fit(X)
Xnew = pca.transform(X)

print(pca.components_.T)
print(np.allclose(Xnew, Y))

# 7. Plot the data
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.title('PCA: eigenvalues')
plt.bar(np.array([0.,0.1]), evals, width=0.1)
plt.xlim(0.,1.)

plt.subplot(1,2,2)
plt.title('PCA: data and eigenvectors')
plt.scatter(X[:,0], X[:,1], color='b')
plt.arrow(0, 0, np.sqrt(evals[0])*P[0,0], np.sqrt(evals[0])*P[1,0],
          length_includes_head = True, head_width = 0.15, color='b')
plt.arrow(0, 0, np.sqrt(evals[1])*P[0,1], np.sqrt(evals[1])*P[1,1],
          length_includes_head = True, head_width = 0.15, color='b')
plt.scatter(Y[:,0], Y[:,1], color='r')
plt.arrow(0, 0, np.sqrt(evals[0]), 0, length_includes_head = True, head_width = 0.15, color='r')
plt.arrow(0, 0, 0, np.sqrt(evals[1]), length_includes_head = True, head_width = 0.15, color='r')
plt.show()

The covariance matrix can be recovered from these eigenvalues and eigenvectors.

In [ ]:

# Using the eigenvectors and eigenvalues, we can of course recover the covariance matrix
# Note: if P is not square, we have to fill it up.
np.allclose(CovX, P.dot(np.diag(evals)).dot(P.T))

Few properties:

• Applying PCA several times is the same as applying it one time. This is because PCA diagonalizes our matrix, thus applying PCA on a diagonal matrix will result in the same matrix.
• Applying PCA on a part of the data and another PCA on the other part, then applying PCA on the concatenation of both do not result in the same matrix as applying PCA on the whole data.
• Applying PCA on a rotated matrix does not give the same result as applying PCA on the initial matrix.

In [ ]:

# Here is the general method

def pca(X, normalize=False, copy=True):
    if copy:
        X = np.copy(X)

    # 1. Center the data
    mean = X.mean(axis=0)
    X -= mean
    N = X.shape[0]
    
    if normalize:
        X /= X.std(axis=0)

    # 2. Compute the covariance matrix
    CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))

    # 3. Compute the eigenvectors of this covariance matrix
    # np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix
    evals, evecs = np.linalg.eigh(CovX)

    # 4. Sort the eigenvalues (in decreasing order) and eigenvectors
    idx = np.argsort(evals)[::-1]
    evals, evecs = evals[idx], evecs[:,idx]

    return evals, evecs

Applying PCA on time series

A fundamental question when applying PCA on time series is how to visualize this high dimensional data. Indeed, a sample \pmb{x}(t) \in \mathbb{R}^T. One way is to plot this \pmb{x}(t) where the x-axis is the time, and y-axis is x(t). Each time t_i (\forall i \in \{0,...,T\}) represents a dimension. By plotting a vertical line at time t=t_i, we can see the variance in this particular dimension.

For an infinite vector, or function, check about Functional PCA.

Recursive PCA

Let's now apply recursive PCA on this toy example, with 3 new samples coming at different time steps.

In [ ]:

# Let's augment our matrix with 3 new samples
Xs = np.array([[-1,1],
               [-3,0],
               [-4,-1]], dtype=np.float64)

n_components = 2
k = float(N)
R = CovX
X_aug = X
mean = X.mean(axis=0).reshape(-1,1)
print(evals)
for x in Xs:
    x = x.reshape(-1,1)
    X_aug = np.vstack((X_aug, x.T))
    X_aug1 = X_aug - X_aug.mean(axis=0)
    pca = PCA(n_components=n_components)
    pca = pca.fit(X_aug1)
    #print(pca.get_covariance())

    # Recursive PCA
    new_mean = k/(k+1) * mean + 1./(k+1) * x
    diff_mean = (new_mean - mean)
    R = (k-1)/k * R + diff_mean.dot(diff_mean.T) + 1./k * (x-new_mean).dot((x-new_mean).T)
    #print(R)
    print("The cov of the whole augmented matrix is equal to the recursive cov: {0}".format(
            np.allclose(pca.get_covariance(), R)))
    k+=1
    mean = new_mean
    
    evals = np.linalg.eigh(R)[0]
    idx = np.argsort(evals)[::-1]
    evals = evals[idx]
    print(evals)

# Compute the new projection matrix
evals, evecs = np.linalg.eigh(R)
idx = np.argsort(evals)[::-1]
evals, evecs = evals[idx], evecs[:,idx]
P = evecs[:,:n_components]
Y = X_aug.dot(P)
    
# Plot the data
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.title('PCA: eigenvalues')
plt.bar(np.array([0.,0.1]), evals, width=0.1)
plt.xlim(0.,1.)

plt.subplot(1,2,2)
plt.title('PCA: data and eigenvectors')
plt.scatter(X_aug[:,0], X_aug[:,1], color='b')
plt.arrow(0, 0, np.sqrt(evals[0])*P[0,0], np.sqrt(evals[0])*P[1,0],
          length_includes_head = True, head_width = 0.15, color='b')
plt.arrow(0, 0, np.sqrt(evals[1])*P[0,1], np.sqrt(evals[1])*P[1,1],
          length_includes_head = True, head_width = 0.15, color='b')
plt.scatter(Y[:,0], Y[:,1], color='r')
plt.arrow(0, 0, np.sqrt(evals[0]), 0, length_includes_head = True, head_width = 0.15, color='r')
plt.arrow(0, 0, 0, np.sqrt(evals[1]), length_includes_head = True, head_width = 0.15, color='r')
plt.show()

Let's now add a sample that can not be modeled by a linear combination of the principal axes, i.e. which is orthogonal to the current covariance matrix. Then, as usual, let's apply recursive PCA on it.

In [ ]:

# Let's add a sample that can not be modeled by a linear combination of the principal axes
# i.e. which is orthogonal to the current covariance matrix.
# Then, let's apply recursive PCA on it.

# New 3D sample
x = np.array([1,-1,1], dtype=np.float64).reshape(-1,1)

# Reshaping previous values (pad a column/row of zeros)
X_aug = np.pad(X_aug, ((0,0), (0,1)), mode='constant') # Nx(T+1)
mean = np.pad(mean, ((0,1),(0,0)), mode='constant')
R = np.pad(R, ((0,1),(0,1)), mode='constant')

# Adding new sample and compute mean
X_aug = np.vstack((X_aug, x.T))
X_aug1 = X_aug - X_aug.mean(axis=0)

# Use sklearn PCA
n_components = 3
pca = PCA(n_components=n_components)
pca = pca.fit(X_aug1)
#print(pca.get_covariance())

# Recursive PCA
new_mean = k/(k+1) * mean + 1./(k+1) * x
diff_mean = (new_mean - mean)
R = (k-1)/k * R + diff_mean.dot(diff_mean.T) + 1./k * (x-new_mean).dot((x-new_mean).T)
#print(R)
print("The cov of the whole augmented matrix is equal to the recursive cov: {0}".format(
            np.allclose(pca.get_covariance(), R)))
k+=1
mean = new_mean
#print('-'*30)

# Compute the new projection matrix
evals, evecs = np.linalg.eigh(R)
idx = np.argsort(evals)[::-1]
evals, evecs = evals[idx], evecs[:,idx]
P = evecs[:,:n_components]
Y = X_aug.dot(P)
    
# Plot the data
fig = plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.title('PCA: eigenvalues')
plt.bar(np.array([0.,0.1,0.2]), evals, width=0.1)
plt.xlim(0.,1.)

ax = fig.add_subplot(122, projection='3d')
ax.set_title('PCA: data and eigenvectors')
ax.scatter(X_aug[:,0], X_aug[:,1], X_aug[:,2])
# From https://stackoverflow.com/questions/22867620/putting-arrowheads-on-vectors-in-matplotlibs-3d-plot
for v in evecs:
    a = Arrow3D([0., v[0]], [0., v[1]], [0., v[2]],
                mutation_scale=20, lw=1, arrowstyle="-|>", color="b")
    ax.add_artist(a)
ax.scatter(Y[:,0], Y[:,1], Y[:,2], color='r')
a = Arrow3D([0., evals[0]], [0., evals[1]], [0., evals[2]],
            mutation_scale=20, lw=1, arrowstyle="-|>", color="r")
ax.add_artist(a)

plt.show()

PCA as an Optimization Problem

PCA can be viewed as an optimization problem in 2 different ways. Theses 2 approaches are equivalent.

1. Maximize the variance of the projected data.
2. Minimize the reconstruction error in a least-square sense.

Mathematically, here is the optimization problem that we are trying to solve:

\begin{equation} \max_{\pmb{v_i}} : \pmb{v_i}^T \pmb{X}^T \pmb{X v_i} \quad \mbox{subj. to} \quad \begin{array}{l} \pmb{v_i}^T \pmb{v_i} = 1 \ \pmb{v_i}^T \pmb{v_j} = 0 \end{array}, \end{equation} \forall i \in \{1,...,D\}, \forall 1\leq j < i.

Nice references:

• What is the objective fct of PCA? (StackExchange)
• PCA objective function: what is the connection between maximizing variance and minimizing error? (StackExchange)
• Everything you did and didn't know about PCA (blog)
• "PCA and Optimization - A Tutorial", 2015 (paper)

In [ ]:

# data
# Toy data
X = np.array([[-1, -1],
              [-2, -1],
              [-3, -2],
              [1, 1],
              [2, 1],
              [3, 2]], dtype=np.float64) # NxT
N = X.shape[0]

# PCA
evals, evecs = pca(X)
print(evals)
print(evecs)

Using Scipy

In [ ]:

from scipy.optimize import minimize

# PCA as an optimization

# cache the 'covariance' matrix
C = X.T.dot(X)/(N-1)

# define objective function to MINIMIZE
f = lambda u: -(u.T.dot(C)).dot(u)

# define initial guess
x0 = np.zeros((2,1))

# define optimization method
# By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds
method = None

# define constraints
constraints = [{'type': 'eq', 'fun': lambda u: u.T.dot(u) - 1}]

# Minimize --> it returns an instance of OptimizeResult
u1 = minimize(f, x0, method=method, constraints=constraints)
#print(u1)
u1 = u1.x.reshape(-1,1) # get solution

# Add constraint
constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})
u2 = minimize(f, x0, method=method, constraints=constraints)
#print(u2)
u2 = u2.x.reshape(-1,1)

# stack the optimized vector found
P = np.hstack((u1,u2))
print(P)

# Plot
plt.title('PCA: data and eigenvectors')
plt.scatter(X[:,0], X[:,1], color='b')
plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')
plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')
plt.show()

In [ ]:

# Define PCA optimization method
def pca_scipy(X, rough_param=0.0, normalize=False, copy=True):
    """
    Compute PCA on the given data using an optimization process.
    """
    if copy:
        X = np.copy(X)

    # center the data
    mean = X.mean(axis=0)
    X -= mean
    N = X.shape[0]
    T = X.shape[1]

    # normalize
    if normalize:
        X /= X.std(axis=0)

    
    # compute 'covariance' matrix and cache it
    N = X.shape[0]
    C = X.T.dot(X)/(N-1)

    # define objective function to MINIMIZE
    #f = lambda u: -(u.T.dot(C)).dot(u)
    def f(u):
        pen = 0
        if rough_param != 0 and u.size > 2:
            ddu = np.diff(np.diff(u))
            rough_pen = (ddu**2).sum()
            pen = rough_param*rough_pen
            #if u.size > 4:
            #    ddddu = np.diff(np.diff(ddu))
            #    rough_pen = (ddddu**2).sum()
            #    pen += rough_param*rough_pen
        loss = -(u.T.dot(C)).dot(u)
        return loss + pen

    # define initial guess
    x0 = np.ones((T,)) #np.zeros((T,))

    # define optimization method
    # By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds
    # If constraints, it will be 'SLSQP' (Sequential Least SQuares Programming)
    # 'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'Newton-CG', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP', 'dogleg', 'trust-ncg'
    # 'Nelder-Mead', 'Powell', 'CG', 'Newton-CG', 'TNC', 'COBYLA', 'dogleg', 'trust-ncg' can not handle (eq) constraints
    # 'BFGS', 'L-BFGS-B' do not work
    method = 'SLSQP'

    # define 1st constraints: norm of 1
    constraints = [{'type': 'eq', 'fun': lambda u: u.T.dot(u) - 1}]

    # optimize recursively
    evals, evecs = [], []
    messages = {}
    for i in range(T):
        if i != 0:
            # add orthogonality constraint
            constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})

        # minimize --> it returns an instance of OptimizeResult
        u1 = minimize(f, x0, method=method, constraints=constraints)
        if not u1.success:
            messages[i] = u1.message

        # get 'eigenvalue'
        evals.append(-u1.fun)

        # get solution
        u1 = u1.x
        evecs.append(u1)

    return np.array(evals), np.array(evecs).T, messages

In [ ]:

evals, evecs, messages = pca_scipy(X)
P = evecs
print(messages)
print(evals)
print(P)

# Plot
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.title('PCA: eigenvalues')
plt.bar(np.array([0.,0.1]), evals, width=0.1)
plt.xlim(0.,1.)

plt.subplot(1,2,2)
plt.title('PCA: data and eigenvectors')
plt.scatter(X[:,0], X[:,1], color='b')
plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')
plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')
plt.show()

Using CVXPY

Note: You cannot use cvxpy for this problem, as we are trying to maximize a convex function, and cvxpy only accepts to maximize a concave fct.

In [ ]:

import cvxpy as cvx

# cache the 'covariance' matrix
C = X.T.dot(X)/(N-1)

# define vector to optimize
u1 = cvx.Variable(X.shape[1])
print(cvx.quad_form(u1, C).is_dcp())
print(cvx.quad_form(u1, C).is_quadratic())

# define objective fct to maximize
#f = cvx.Maximize(u1.T*C*u1) 
f = cvx.Maximize(cvx.quad_form(u1, C)) # this does not work!
#f = cvx.Minimize(cvx.quad_form(u1, C)) # this works (if no constraints) but that is not what we want to achieve!
constraints = [u1.T*u1 == 1]
prob = cvx.Problem(f, constraints)

result = prob.solve()
print(u1.value)

Using NLopt

Nonlinear optimization algorithms that can handle nonlinear inequality and equality constraints are:

• ISRES (Improved Stochastic Ranking Evolution Strategy) \rightarrow global derivative-free
• COBYLA (Constrained Optimization BY Linear Approximations) \rightarrow local derivative-free
• SLSQP (Sequential Least-SQuares Programming) \rightarrow local gradient-based
• AUGLAG (AUGmented LAGrangian) \rightarrow global/local derivative-free/gradient based (determined based on the subsidiary algo)

More information about the various algorithms can be found on this link.

In [ ]:

# Playground with NLopt
import nlopt

nlopt_results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',
                 5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',
                 -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}
n_iter = 0
N = X.shape[0]

# cache the 'covariance' matrix
C = X.T.dot(X)/(N-1)

# define random seed
nlopt.srand(125)

# define which solver to use
#optimizer = nlopt.GN_ISRES
#optimizer = nlopt.LN_COBYLA
optimizer = nlopt.LD_SLSQP
#optimizer = nlopt.LD_AUGLAG # nlopt.AUGLAG, nlopt.AUGLAG_EQ, nlopt.LD_AUGLAG,
                            # nlopt.LD_AUGLAG_EQ, nlopt.LN_AUGLAG, nlopt.LN_AUGLAG_EQ

# if nlopt.AUGLAG, we have to define a subsidiary algo
suboptimizer = nlopt.LD_SLSQP #nlopt.LN_COBYLA

# define objective function
def f(x, grad):
    global n_iter
    n_iter += 1
    if grad.size > 0:
        grad[:] = 2*x.T.dot(C)
    return x.T.dot(C).dot(x)

# define norm constraint
def c1(x, grad):
    if grad.size > 0:
        grad[:] = 2*x
    return (x.T.dot(x) - 1)

# create optimizer
n = X.shape[1] # nb of parameters to optimize, size of the vector
opt = nlopt.opt(optimizer, n)
print("Algo: %s" % opt.get_algorithm_name())
opt.set_max_objective(f)

# if nlopt.GN_ISRES, we can define the population size
opt.set_population(0) # by default for ISRES: pop=20*(n+1)

# if nlopt.AUGLAG, set the subsidiary algo
subopt = nlopt.opt(suboptimizer, n)
subopt.set_lower_bounds(-1)
subopt.set_upper_bounds(1)
#subopt.set_ftol_rel(1e-2)
#subopt.set_maxeval(100)
opt.set_local_optimizer(subopt)

# define bound constraints (should be between -1 and 1 because the norm should be 1)
opt.set_lower_bounds(-1.)
opt.set_upper_bounds(1.)

# define equality constraints
opt.add_equality_constraint(c1, 0)
#opt.add_equality_mconstraint(constraints, tol)

# define stopping criteria
#opt.set_stopval(stopval)
opt.set_ftol_rel(1e-8)
#opt.set_xtol_rel(1e-4)
opt.set_maxeval(100000) # nb of iteration
opt.set_maxtime(10) # time in secs

# define initial value
#x0 = np.zeros((n,))
x0 = np.array([0.1,0.1])

# optimize
x = x0
try:
    x = opt.optimize(x0)
except nlopt.RoundoffLimited as e:
    pass
max_value = opt.last_optimum_value()
result = opt.last_optimize_result()

print("Max value: %f" % max_value)
print("Nb of iterations: %d" % n_iter)
print("Result: %s" % nlopt_results[result])
print("Optimized array:")
print(x)

In [ ]:

# Define PCA optimization method formally using NLopt
def center_data(X, normalize=False, copy=True):
    if copy:
        X = np.copy(X)

    # center the data
    mean = X.mean(axis=0)
    X -= mean
    N = X.shape[0]
    T = X.shape[1]

    # normalize
    if normalize:
        X /= X.std(axis=0)
    
    return X

class OrthogonalConstraint(object):
    
    def __init__(self, v):
        self.v = np.copy(v)
        
    def constraint(self, x, grad):
        if grad.size > 0:
            grad[:] = self.v
        return (x.T.dot(self.v))
        

def pca_nlopt(X, method=None, submethod=None, rough_param=0.0, normalize=False, copy=True):
    """
    Compute PCA on the given data using nlopt.
    
    :param (str) method: it can take the following value: 'SLSQP', 'ISRES',
                         'COBYLA', 'AUGLAG'. By default, it will be 'SLSQP'.
    :param (str) submethod: this needs to be defined only if method is 'AUGLAG'.
                            By default, it will be 'SLSQP'.
    """
    # center the data
    X = center_data(X, normalize=normalize, copy=copy)

    # compute 'covariance' matrix and cache it
    N = X.shape[0]
    C = X.T.dot(X) / (N-1)
    
    # define useful variables
    nlopt_results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',
                 5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',
                 -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}
    n = X.shape[1] # nb of parameters to optimize, size of the vector
    
    # define random seed
    nlopt.srand(125)

    # define which solver (and possibly subsolver) to use
    def get_opt(method):
        if method == 'ISRES':
            return nlopt.opt(nlopt.GN_ISRES, n)
        elif method == 'COBYLA':
            return nlopt.opt(nlopt.LN_COBYLA, n)
        elif method == 'SLSQP':
            return nlopt.opt(nlopt.LD_SLSQP, n)
        elif method == 'AUGLAG':
            return nlopt.opt(nlopt.AUGLAG, n)
        else:
            raise NotImplementedError("The given method has not been implemented")

    if method is None:
        method = 'SLSQP'            
    opt = get_opt(method)
    if method == 'AUGLAG':
        if submethod is None:
            submethod = 'SLSQP'
        elif submethod == 'AUGLAG':
            raise ValueError("Submethod should be different from AUGLAG")
        subopt = get_opt(submethod)
        subopt.set_lower_bounds(-1)
        subopt.set_upper_bounds(1)
        #subopt.set_ftol_rel(1e-2)
        #subopt.set_maxeval(100)
        opt.set_local_optimizer(subopt)
        
    # define objective function
    def f(x, grad):
        if grad.size > 0:
            grad[:] = 2*x.T.dot(C)
        return x.T.dot(C).dot(x)

    # define norm constraint
    def c1(x, grad):
        if grad.size > 0:
            grad[:] = 2*x
        return (x.T.dot(x) - 1)
    
    # define objective function
    opt.set_max_objective(f)
        
    # if nlopt.GN_ISRES, we can define the population size
    opt.set_population(0) # by default for ISRES: pop=20*(n+1)
    
    # define bound constraints (should be between -1 and 1 because the norm should be 1)
    opt.set_lower_bounds(-1.)
    opt.set_upper_bounds(1.)

    # define equality constraints
    opt.add_equality_constraint(c1, 0)
    #opt.add_equality_mconstraint(constraints, tol)

    # define stopping criteria
    #opt.set_stopval(stopval)
    opt.set_ftol_rel(1e-8)
    #opt.set_xtol_rel(1e-4)
    opt.set_maxeval(100000) # nb of iteration
    opt.set_maxtime(2) # time in secs

    # define initial value
    x0 = np.array([0.1]*n) # important that the initial value ≠ 0 for the computation of the grad!

    evals, evecs, msgs = [], [], {}
    for i in range(n):
        if i > 0:
            c = OrthogonalConstraint(x)
            opt.add_equality_constraint(c.constraint, 0)
        # optimize
        try:
            x = opt.optimize(x0)
        except nlopt.RoundoffLimited as e:
            pass

        # save values
        evecs.append(x)
        evals.append(opt.last_optimum_value())
        msgs[i] = nlopt_results[opt.last_optimize_result()]

    return np.array(evals), np.array(evecs).T, msgs

In [ ]:

method = 'SLSQP' # 'SLSQP', 'COBYLA', 'ISRES', 'AUGLAG'
submethod = None

evals, P, msgs = pca_nlopt(X, method=method, submethod=submethod)
print(msgs)
print(evals)
print(P)

# Plot
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.title('PCA: eigenvalues')
plt.bar(np.array([0.,0.1]), evals, width=0.1)
plt.xlim(0.,1.)

plt.subplot(1,2,2)
plt.title('PCA: data and eigenvectors')
plt.scatter(X[:,0], X[:,1], color='b')
plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')
plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')
plt.show()

For comparison purpose, we obtained the following values with scipy ('SLSQP'):

[ 7.93954312 0.06045688]
[[ 0.83849224 -0.54491355]
[ 0.54491353 0.83849226]]

And these values using std PCA:

[ 7.93954312 0.06045688]
[[-0.83849224 0.54491354]
[-0.54491354 -0.83849224]]

Using IPopt

In [ ]:

# Playground with IPopt
import ipopt

# define useful vars
n = X.shape[1]
N = X.shape[0]
C = X.T.dot(X) / (N-1)

# define initial value
x0 = np.array([0.1]*n)

# define (lower and upper) bound constraints
lb = [-1]*n
ub = [1]*n

# define constraints
cl = [1]
cu = [1]

class Opt(object):
    
    def __init__(self, verbose=True):
        self.verbose = verbose
        self.iter_count = 0

    def objective(self, x):
        # objective fct to minimize
        return -x.T.dot(C).dot(x)
    
    def gradient(self, x):
        # grad of the objective fct
        return -2*x.T.dot(C)
    
    def constraints(self, x):
        # norm constraint
        return x.T.dot(x)
    
    def jacobian(self, x):
        return 2*x
    
    #def hessian(self, x):
    #    pass
    
    def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,
                     regularization_size, alpha_du, alpha_pr, ls_trials):
        if self.verbose:
            print("Objective value at iteration #%d: %g" % (iter_count, obj_value))
        self.iter_count = iter_count

opt = Opt(verbose=False)
nlp = ipopt.problem(n=n, m=len(cl), problem_obj=opt, lb=lb, ub=ub, cl=cl, cu=cu)

x, info = nlp.solve(x0)
print("Max value: %f" % -info['obj_val'])
print("Nb of iterations: %d" % opt.iter_count)
print("Result: %s" % info['status_msg'])
print("Optimized array:")
print(x)

In [ ]:

# Define PCA optimization method formally using IPopt
class NormConstraint(object):
    
    def __init__(self):
        pass
    
    def constraint(self, x):
        return x.T.dot(x)
    
    def jacobian(self, x):
        return 2*x
    
class OrthogonalConstraint(object):
    
    def __init__(self, v):
        self.v = np.copy(v)
    
    def constraint(self, x):
        return x.T.dot(self.v)
    
    def jacobian(self, x):
        return self.v

class IPopt(object):

    def __init__(self, verbose=True):
        self.verbose = verbose
        self.iter_count = 0
        self.cons = []
    
    def add_constraint(self, c):
        self.cons.append(c)

    def objective(self, x):
        # objective fct to minimize
        return -x.T.dot(C).dot(x)
    
    def gradient(self, x):
        # grad of the objective fct
        return -2*x.T.dot(C)
    
    def constraints(self, x):
        return np.array([c.constraint(x) for c in self.cons])
    
    def jacobian(self, x):
        return np.array([c.jacobian(x) for c in self.cons])
    
    #def hessian(self, x):
    #    pass
    
    def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,
                     regularization_size, alpha_du, alpha_pr, ls_trials):
        if self.verbose:
            print("Objective value at iteration #%d: %g" % (iter_count, obj_value))
        self.iter_count = iter_count

def pca_ipopt(X, rough_param=0.0, normalize=False, copy=True):
    """
    Compute PCA on the given data using ipopt.
    """
    # center the data
    X = center_data(X, normalize=normalize, copy=copy)
    
    # define useful vars
    n = X.shape[1]
    N = X.shape[0]
    C = X.T.dot(X) / (N-1)

    # define initial value
    x0 = np.array([0.1]*n)

    # define (lower and upper) bound constraints
    lb = [-1]*n
    ub = [1]*n

    # define constraints
    cl = [1] + [0]*(n-1)
    cu = [1] + [0]*(n-1)

    evals, evecs, msgs = [], [], {}
    opt = IPopt(verbose=False)
    opt.add_constraint(NormConstraint())
    for i in range(n):
        if i > 0:
            opt.add_constraint(OrthogonalConstraint(x))
        
        i1 = i+1
        nlp = ipopt.problem(n=n, m=len(cl[:i1]), problem_obj=opt,
                            lb=lb, ub=ub, cl=cl[:i1], cu=cu[:i1])
        
        # solve problem
        x, info = nlp.solve(x0)
        
        evecs.append(x)
        evals.append(-info['obj_val'])
        msgs[i] = info['status_msg']
        
    return np.array(evals), np.array(evecs).T, msgs

In [ ]:

evals, P, msgs = pca_ipopt(X)
print(msgs)
print(evals)
print(P)

# Plot
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.title('PCA: eigenvalues')
plt.bar(np.array([0.,0.1]), evals, width=0.1)
plt.xlim(0.,1.)

plt.subplot(1,2,2)
plt.title('PCA: data and eigenvectors')
plt.scatter(X[:,0], X[:,1], color='b')
plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')
plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')
plt.show()