I am using singular vector decomposition on a matrix and obtaining the U, S and Vt matrices. At this point, I am trying to choose a threshold for the number of dimensions to retain. I was suggested to look at a scree plot but am wondering how to go about plotting it in numpy. Currently, I am doing the following using numpy and scipy libraries in python:

U, S, Vt = svd(A)

Any suggestions?


closed as off-topic by amoeba, Peter Flom Jan 28 at 12:07

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – amoeba, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    $\begingroup$ take the diagonal of S, if it is not already a diagonal, square it, sort it in decreasing order, take the cumulative sum, divide by the last value, then plot it. $\endgroup$ – shabbychef Jul 9 '11 at 4:39
  • $\begingroup$ @shabbychef: You mean, take the cumulative sum and divide by the sum of all the values right? $\endgroup$ – Legend Jul 10 '11 at 1:24
  • $\begingroup$ yes. In matlab, it would be [U,S,V] = svd(X);S = cumsum(sort(diag(S).^2,1,'descend'));S = S ./ S(end);plot(S); $\endgroup$ – shabbychef Jul 10 '11 at 2:35

Here is an example that can be pasted to an IPython prompt and generate an image like below (it uses random data):

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

#Make a random array and then make it positive-definite
num_vars = 6
num_obs = 9
A = np.random.randn(num_obs, num_vars)
A = np.asmatrix(A.T) * np.asmatrix(A)
U, S, V = np.linalg.svd(A) 
eigvals = S**2 / np.sum(S**2)  # NOTE (@amoeba): These are not PCA eigenvalues. 
                               # This question is about SVD.

fig = plt.figure(figsize=(8,5))
sing_vals = np.arange(num_vars) + 1
plt.plot(sing_vals, eigvals, 'ro-', linewidth=2)
plt.title('Scree Plot')
plt.xlabel('Principal Component')
#I don't like the default legend so I typically make mine like below, e.g.
#with smaller fonts and a bit transparent so I do not cover up data, and make
#it moveable by the viewer in case upper-right is a bad place for it 
leg = plt.legend(['Eigenvalues from SVD'], loc='best', borderpad=0.3, 
                 shadow=False, prop=matplotlib.font_manager.FontProperties(size='small'),

enter image description here

  • $\begingroup$ Hermann: +1 Thank You for your time! I know it's been a long time but nevertheless this is really good to have :) $\endgroup$ – Legend Oct 18 '11 at 19:08
  • $\begingroup$ what is num_vars? it doesn't seem to be defined in your script. $\endgroup$ – TheChymera Jul 30 '14 at 17:00
  • $\begingroup$ @TheChymera - Thanks for catching this, I have updated my response. $\endgroup$ – Josh Hemann Aug 12 '14 at 19:01
  • $\begingroup$ @Josh Hemann yeah, i also figured this out in the mean time - but I think it might be better to compute it from the shape of A $\endgroup$ – TheChymera Aug 12 '14 at 19:41
  • 1
    $\begingroup$ @JoshHemann can you explain this: eigvals = S** 2 / np.cumsum(S)[-1] ?? I have seen based on some papers eigvals = S**2 / (n-1) where n is the number of features $\endgroup$ – makis Mar 6 '18 at 19:09

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