# Implementation of a new algorithm for sklearn

In the Python library, sklearn is implemented the algorithm for SparsePCA.

I have written the code for a another version of this algorithm that is much faster in some situations. I have not enough experience with Python and sklearn in order to write code that could be loaded in the sklearn repositories.

It would be nice if you can help me to write this code in order to make it 100% compatible with the sklearn standards. If you have experience with machine learning would be nice if you help me to consider bad situations that I may have not considered in order to make the algorithm more robust.

import numpy as np
import sklearn as sk
import pandas as pd
from sklearn.decomposition import SparsePCA, PCA

import random

class IterativeSPCA():

def __init__(self, Npc=None, alpha=1, max_iter=5000, tol=1e-8):
self.Npc = Npc #number of components
self.alpha = alpha  # penalty value
self.max_iter = max_iter # max number of iteration
self.tol = tol # accepted tollerance error

def spca_iterate(self, u_old, v_old):
# this is the equivalent of the nipals algorithm but designed in order to
v_new = np.zeros(shape=v_old.shape)
u_new = np.zeros(shape=u_old.shape)

for iteration in range(self.max_iter): #repeat the iteration non more than max iter
y = np.dot(self.X.T, u_old).squeeze()
v_new = (np.sign(y)*np.maximum(np.abs(y)-self.alpha, 0)).squeeze()
norm_v = np.linalg.norm(v_new)**2
if norm_v == 0: #if norm v is 0 it means that no components have been selected. So you should reduce the penalty lambda

break

x_times_v = np.dot(self.X, v_new)
x_times_v.shape = (self.m, 1)
u_new = x_times_v/np.linalg.norm(x_times_v)

if np.linalg.norm(v_new)==0:#check again if the norm is 0
break
if np.linalg.norm(v_new - v_old)<self.tol or np.linalg.norm(-v_new - v_old) < self.tol:#check if there is convergence
break
u_old = u_new
v_old = v_new

if iteration == self.max_iter-1:
print("No Convergence. Error!!!")

norm_v = np.linalg.norm(v_new)
v_new.shape = (self.n, 1)
v_final = sk.preprocessing.normalize(v_new, axis=0)
u_final = u_new * norm_v
return v_final, u_final

def fit(self, X):

self.m, self.n = X.shape
self.components_ = np.zeros(shape=(self.n,self.Npc))
self.sT = np.zeros(shape=(self.m,self.Npc))
self.X = X

for i in range(self.Npc):

pca = PCA(n_components=1).fit(self.X)
v = pca.components_
v = sk.preprocessing.normalize(v)
u = self.X.dot(v.T)
s = np.linalg.norm(u)

u_old = u
v_old = v*s

v_final, u_final = self.spca_iterate(u_old=u_old, v_old=v_old)
self.sT[:, i] = u_final.squeeze()
self.components_[:, i] = v_final.squeeze()
u_final.shape = (self.m, 1)
self.X = self.X-np.dot(u_final, v_final.T)
return self.components_, self.sT

print("start")
X = np.random.random((100, 20)) #generate a random matrix
X = sk.preprocessing.scale(X)   # scale the data
alpha = 1 # set the penalty
nPCs = 2 # set the number of components

import time

start_time = time.time()
a = IterativeSPCA(Npc=nPCs, alpha=alpha)
sP, sT = a.fit(X)
time1 = time.time()

spca = SparsePCA(n_components=nPCs, alpha=alpha, ridge_alpha=0)
spca.fit(X)
time2 = time.time()

import pandas as pd
df = pd.DataFrame(columns=['it_spc1','spc1'])
df['it_spc1'] = a.components_[:,0]
df['spc1'] = spca.components_.T[:,0]
df['it_spc2'] = a.components_[:,1]
df['spc2'] = spca.components_.T[:,1]
print("time iterative SPCA=", time1-start_time, "time SparsePCA=", time2-time1)
print("finish!")

• Do you want to contribute this back? I imagine the best thing to do is ask about it on their mailing list or issues page, and ask for a code review. – Veedrac May 26 '14 at 10:04
• Thanks for answering. I have asked on the mailing list but I didn't receive any feedback – Donbeo May 26 '14 at 15:43

Thank you for posting this.

You didn't include any result timings. CPU performance and error performance on an available dataset would be of interest, for evaluating against sklearn's baseline code. Highly correlated inputs would be more interesting than uniform random numbers.

Noticing the import just above it makes the class name perfectly clear, but nonetheless please give it a docstring explaining:

"""Iterative Sparse PCA (principal component analysis)."""


    self.Npc = Npc  # number of principal components


or, better, put it in a docstring. Consider deleting the redundant max_iter comment, or perhaps rephrase it to describe max number of principal components if that is its function.

Typo: tolerance

def spca_iterate(self, u_old, v_old):


The public API should probably talk about u_init & v_init.

Typo: no more (but consider deleting this redundant comment)

    v_new = np.zeros(shape=v_old.shape)
u_new = np.zeros(shape=u_old.shape)


This pair of assignments appears to have no effect.

        if norm_v == 0:  # if norm v is 0 it means that no components have been selected. So you should reduce the penalty lambda
break


This seems to be speaking about this author's alpha in terms of another author's lambda. Your docstring should definitely mention the effect of alpha on the solution's sparsity.

        if np.linalg.norm(v_new)==0:#check again if the norm is 0
break
if np.linalg.norm(v_new - v_old)<self.tol or np.linalg.norm(-v_new - v_old) < self.tol:#check if there is convergence
break


The first test is pointless (redundant). Please delete it. The 2nd test is phrased oddly. Please express it as:

        if abs(np.linalg.norm(v_new - v_old)) < self.tol:  # check for convergence
break


This assignment:

        x_times_v.shape = (self.m, 1)


and the v_new .squeeze(), make me concerned that we didn't take an opportunity to fix up shapes before the loop began. The squeeze returns either a view or the underlying object. A single view would be just fine, but a chain of max_iter views would be a problem. Can we maybe banish the shape concerns to the init code at top of function, outside the loop?

        print("No Convergence. Error!!!")


A raise might be more appropriate here.

    norm_v = np.linalg.norm(v_new)


We were squaring before. This statement might be fine, but it warrants a comment, as one norm_v is not like the other. Or, put another way, the whole method really really needs a specific literature reference describing the method.

It's not actually clear what the squaring was about. If FP negative zero was a concern, then abs() should suffice.

fit() would like a docstring, please.

    self.components_ = np.zeros(shape=(self.n,self.Npc))
...
self.components_[:, i] = v_final.squeeze()


This made no sense to me. It is conventional to use underscore suffix to avoid conflict with pre-existing symbols (e.g. min_) and prefix for private variables. But here it appears you just wanted a local variable (not an attribute) named components. For that matter, it's not clear what's going on with m, n, sT, & X. It appears they should be locals.

X = np.random.random((100, 20)) #generate a random matrix
X = sk.preprocessing.scale(X)   # scale the data


alpha = 1  # set the penalty


Good. Very clear, consistent with the literature.

nPCs = 2  # set the number of components


You felt a need to write a comment. This suggests nPCs was misnamed, and would be clearer as num_components. Similarly for Npc in the ctor.

Consider using a private _spca_iterate name.

The a name helps nobody. Consider turning

a = IterativeSPCA(Npc=nPCs, alpha=alpha)
sP, sT = a.fit(X)


into

sP, sT = IterativeSPCA(Npc=nPCs, alpha=alpha).fit(X)


Maybe a single fit() function, with no class, would suffice.

The sP and sT names could be more helpful. I understand that P denotes principle components, but you didn't offer a literature reference so T is a bit obscure, it doesn't seem to relate to "feature".

df['it_spc1'] = a.components_[:,0]


What is going on here? Didn't you just put that in sP, but now you're unwilling to use it?

df['spc1'] = spca.components_.T[:,0]


Oh, I see. You were going for API compatibility. I guess you know already how I feel about SPCA fit() being void. I cannot fathom why you return transpose of what SPCA offers.

Ok, maybe .components_ is an important part of your public API, at least for backward compatibility. But it seems like you could be friendlier to callers, so it is easier to consume your API and correctly interpret results.