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 
        # obtain a sparse loading
        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


            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
            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
            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

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)
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)
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Veedrac May 26 '14 at 10:04
  • \$\begingroup\$ Thanks for answering. I have asked on the mailing list but I didn't receive any feedback \$\endgroup\$ – 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)."""

Please fix a comment so it reads:

    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

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
        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

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

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

Please delete both remarks.

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)


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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.