Hello Dave,
I don't know if this will be useful to your research, but it may be
worth pointing out in general. As you know PCA (and perhaps some
other spectral algorithms?) use eigenvalues of matrices that can be
factored out as A'A (where ' means transpose). For example, in the
PCA case, if A is a centered data matrix (i.e. the mean value of each
data point has been subtracted off), and if the data elements are row-
vectors, then A'A is the covariance matrix. PCA then examines the
(nonzero) eigenvectors of this covariance matrix A'A.
Now, it's possible to determine the eigenvalues/eigenvectors for A'A
from the matrix AA', which in many useful cases is much smaller than
A'A. For example, imagine that you have 100 data points, each in
10,000 dimensions. (This is common in imaging applications.) A'A is
10,000x10,000, but AA' is only 100x100. We can get the eigenvalues/
vectors of A'A from those of AA', as described below. (This works in
part because in such cases, the larger matrix is rank-deficient, so
there will only be e.g. 100 nonzero eigenvalues anyway.)
If you're lucky enough to have this kind of structure in your
problem, feel free to use the following code which exploits that (and
explains in a bit more detail how it works).
Zach Pincus
Program in Biomedical Informatics and Department of Biochemistry
Stanford University School of Medicine
def _symm_eig(a):
"""Return the eigenvectors and eigenvalues of the symmetric matrix
a'*a. If
a has more columns than rows, then that matrix will be rank-
deficient,
and the non-zero eigenvalues and eigenvectors can be more easily
extracted
from the matrix a*a', from the properties of the SVD:
if a of shape (m,n) has SVD u*s*v', then:
a'*a = v*s'*s*v'
a*a' = u*s*s'*u'
That is, v contains the eigenvectors of a'*a, with s'*s the
eigenvalues,
according to the eigen-decomposition theorem.
Now, let s_hat, an array of shape (m,n), be such that s * s_hat = I
(m,m)
and s_hat * s = I(n,n). With that, we can solve for u or v in
terms of the
other:
v = a'*u*s_hat'
u = a*v*s_hat
"""
m, n = a.shape
if m >= n:
# just return the eigenvalues and eigenvectors of a'a
vecs, vals = _eigh(numpy.dot(a.transpose(), a))
vecs = numpy.where(vecs < 0, 0, vecs)
return vecs, vals
else:
# figure out the eigenvalues and vectors based on aa', which is
smaller
sst_diag, u = _eigh(numpy.dot(a, a.transpose()))
# in case due to numerical instabilities we have sst_diag < 0
anywhere,
# peg them to zero
sst_diag = numpy.where(sst_diag < 0, 0, sst_diag)
# now get the inverse square root of the diagonal, which will
form the
# main diagonal of s_hat
err = numpy.seterr(divide='ignore', invalid='ignore')
s_hat_diag = 1/numpy.sqrt(sst_diag)
numpy.seterr(**err)
s_hat_diag = numpy.where(numpy.isfinite(s_hat_diag), s_hat_diag, 0)
# s_hat_diag is a list of length m, a'u is (n,m), so we can just
use
# numpy's broadcasting instead of matrix multiplication, and
only create
# the upper mxm block of a'u, since that's all we'll use anyway...
v = numpy.dot(a.transpose(), u[:,:m]) * s_hat_diag
return sst_diag, v
def _eigh(m):
"""Return eigenvalues and corresponding eigenvectors of the hermetian
array m, ordered by eigenvalues in decreasing order.
Note that the numpy.linalg.eigh makes no order guarantees."""
values, vectors = numpy.linalg.eigh(m)
order = numpy.flipud(values.argsort())
return values[order], vectors[:,order]
On May 13, 2007, at 8:35 PM, Dave P. Novakovic wrote:
> There are definitely elements of spectral graph theory in my research
> too. I'll summarise
>
> We are interested in seeing the each eigenvector from svd can
> represent in a semantic space
> In addition to this we'll be testing it against some algorithms like
> concept indexing (uses a bipartitional k-meansish method for dim
> reduction)
> also testing against Orthogonal Locality Preserving indexing, which
> uses the laplacian of a similarity matrix to calculate projections of
> a document (or term) into a manifold.
>
> These methods have been implemented and tested for document
> classification, I'm interested in seeing their applicability to
> modelling semantics with a system known as Hyperspace to analog
> language.
>
> I was hoping to do svd to my HAL built out of reuters, but that was
> way too big. instead i'm trying with the traces idea i mentioned
> before (ie contextually grepping a keyword out of the docs to build a
> space around it.)
>
> Cheers
>
> Dave
>
> On 5/14/07, Charles R Harris <[EMAIL PROTECTED]> wrote:
>>
>>
>> On 5/13/07, Dave P. Novakovic <[EMAIL PROTECTED]> wrote:
>>>> Are you trying some sort of principal components analysis?
>>>
>>> PCA is indeed one part of the research I'm doing.
>>
>> I had the impression you were trying to build a linear space in
>> which to
>> embed a model, like atmospheric folk do when they try to invert
>> spectra to
>> obtain thermal profiles. Model based compression would be another
>> aspect of
>> this. I wonder if there aren't some algorithms out there for this
>> sort of
>> thing.
>>
>> Chuck
>>
>>
>> _______________________________________________
>> Numpy-discussion mailing list
>> Numpy-discussion@scipy.org
>> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
