Thread Subject: Computing ( X' * X ) Faster?

I want to compute a real matrix transpose times itself. Is there a faster way to do this than simply computing XtX = X' * X ; ? My X matrix is rather large and not square ( M x N with M > N and M and N are in the tens or even hundreds of thousands).

I would think there might be some shortcuts because XtX will be diagonal, so perhaps there is a way to cut the time approximately in half? I do believe my matrix is also sparse if that might help.

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdupu7$aa6$1@fred.mathworks.com>...
> Is there a faster way to do this than simply computing XtX = X' * X ; ?

I don't think so.

Bruno

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdupu7$aa6$1@fred.mathworks.com>...
> I want to compute a real matrix transpose times itself. Is there a faster way to do this than simply computing XtX = X' * X ; ? My X matrix is rather large and not square ( M x N with M > N and M and N are in the tens or even hundreds of thousands).
>
> I would think there might be some shortcuts because XtX will be diagonal

Is there something missing in your problem statement? Why do you think X' * X will be diagonal?

James Tursa

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdupu7$aa6$1@fred.mathworks.com>...

> I would think there might be some shortcuts because XtX will be diagonal,
====

If this is indeed the case (like James I wonder why), then you can compute the diagonal of XtX according to

d=sum(X.^2)

Since d is the only non-trivial data, you're basically done. If you actually need d loaded into a sparse matrix, then do

A=spdiags(d,0,speye(numel(d)));

> Is there something missing in your problem statement? Why do you think X' * X will >be diagonal?

sorry, i meant to say "symmetric", so that XtX(i,j) = XtX(j,i).

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hduv38$87p$1@fred.mathworks.com>...
>
> > Is there something missing in your problem statement? Why do you think X' * X will >be diagonal?
>
> sorry, i meant to say "symmetric", so that XtX(i,j) = XtX(j,i).

Then there is probably no simplification, barring any special structure in X that you haven't told us about.

It might be worth telling us what you want to do with XtX. It's very unusual to need to compute such a matrix explicitly. For example, if all you want is some matrix-vector product y=XtX*b, then you would certainly compute this according to

y= ( (X*b)' *X)';

without ever having to compute XtX or even to compute Xt.

>
> It might be worth telling us what you want to do with XtX. It's very unusual to need to compute such a matrix explicitly. For example, if all you want is some matrix-vector product y=XtX*b, then you would certainly compute this according to
>

I often find myself computing X'*X, followed by eigenvectors via eigs(XtX), to do PCA. Is there a way to accomplish PCA without X'*X? I've profiled my code, and X'*X is >>50% of my execution time. I'm running under 32-bit Vista, so princomp() or svds() blow up with my 16+k x 2k data matrices.
Thanks.

> It might be worth telling us what you want to do with XtX. It's very unusual to >need to compute such a matrix explicitly. For example, if all you want is some matrix-?vector product y=XtX*b, then you would certainly compute this according to
>
> y= ( (X*b)' *X)';
>
> without ever having to compute XtX or even to compute Xt.

Active set strategies for least squares with non-negativity constraints usually require this, and they appear to be superior in terms of speed to lsqnonneg() that ships with MATLAB. For example, take a look at the NNLS code of Rasmus Bro on the file exchange...

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdv24o$i7d$1@fred.mathworks.com>...
> > It might be worth telling us what you want to do with XtX. It's very unusual to >need to compute such a matrix explicitly. For example, if all you want is some matrix-?vector product y=XtX*b, then you would certainly compute this according to
> >
> > y= ( (X*b)' *X)';
> >
> > without ever having to compute XtX or even to compute Xt.
>
> Active set strategies for least squares with non-negativity constraints usually require this, and they appear to be superior in terms of speed to lsqnonneg() that ships with MATLAB. For example, take a look at the NNLS code of Rasmus Bro on the file exchange...

Note that:

Active set with *box constraint* (such as non-negative) requires only matrix-vector product.

Active set without box constraints (general linear constraints) can be transformed into box-constraint using dual variable. But usually the sparsity of lost because the inverse of the Hessian is needed.

Bruno

"Chad " <parishcm@ornl.gov> wrote in message <hdv1lm$iqp$1@fred.mathworks.com>...
>
> >
> > It might be worth telling us what you want to do with XtX. It's very unusual to need to compute such a matrix explicitly. For example, if all you want is some matrix-vector product y=XtX*b, then you would certainly compute this according to
> >
>
> I often find myself computing X'*X, followed by eigenvectors via eigs(XtX), to do PCA. Is there a way to accomplish PCA without X'*X? I've profiled my code, and X'*X is >>50% of my execution time. I'm running under 32-bit Vista, so princomp() or svds() blow up with my 16+k x 2k data matrices.
> Thanks.

svds would have been my suggestion, but I'm surprised that resorting to X'*X doesn't give you numerical problems...

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hduv38$87p$1@fred.mathworks.com>...
>
> > Is there something missing in your problem statement? Why do you think X' * X will >be diagonal?
>
> sorry, i meant to say "symmetric", so that XtX(i,j) = XtX(j,i).

For full matrices, MATLAB uses BLAS library calls to do matrix multiplies. In the case of a general matrix multiply X' * Y where X and Y are different, the routine DGEMM is called. But if MATLAB recognizes a symmetric form of the multiply, as in your case with X' * X where the two operands are the same, MATLAB instead calls the routines DSYRK and/or DSYR2K that work specifically with this form of multiply, only doing the lower or upper half. Then MATLAB fills in the other half with a copy. This cuts the computation time roughly in 1/2 for these cases. So you are already getting the benefit of the symmetric BLAS routines when you do the X' * X multiply and I don't think you can do anything else to speed it up.

I am not familiar with the sparse matrix multiply routines, but from simple timing tests it doesn't look like there is any special symmetric algorithm involved. e.g.,

>> A = rand(2000,2000);
>> B = rand(2000,2000);
>> tic;A'*B;toc
Elapsed time is 1.139315 seconds.
>> tic;A'*A;toc
Elapsed time is 0.683641 seconds.
>> A = sprand(5000,5000,.01);
>> B = sprand(5000,5000,.01);
>> tic;A'*B;toc
Elapsed time is 0.628661 seconds.
>> tic;A'*A;toc
Elapsed time is 0.636321 seconds.

James Tursa

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hdv2t7$6gr$1@fred.mathworks.com>...

> Active set with *box constraint* (such as non-negative) requires only matrix-vector product.
===============

This is because you can just use the active constraints to eliminate the constrained
x(i) from the quadratic cost function?

  

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdv24o$i7d$1@fred.mathworks.com>...
> > It might be worth telling us what you want to do with XtX. It's very unusual to >need to compute such a matrix explicitly. For example, if all you want is some matrix-?vector product y=XtX*b, then you would certainly compute this according to
> >
> > y= ( (X*b)' *X)';
> >
> > without ever having to compute XtX or even to compute Xt.
>
> Active set strategies for least squares with non-negativity constraints usually require this, and they appear to be superior in terms of speed to lsqnonneg() that ships with MATLAB. For example, take a look at the NNLS code of Rasmus Bro on the file exchange...

=====================

Just to throw another possible approach into the mix, the code below is a Majorize-Minimize algorithm for minimizing

norm(A*x-b).^2 subject to x>=0

It has the advantage that you never have to compute A'*A explicitly, but its performance is normally better the more concentrated AtA is along the diagonal.


W=sum(A.^2,2); %Hessian diagonal
W=1./W; %invert



 Curvs=A'*(sum(A,2).*W); %majorizing curvatures

 x=x0; %Initialize x with some vector x0>=0

for i=1:NumIterations %MAIN LOOP

 gradient=A.'*((Ax-b).*W);
 x=x-gradient./Curvs;
 x(x<0)=0;

end

"Matt " <xys@whatever.com> wrote in message <hdv5d5$d1a$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hdv2t7$6gr$1@fred.mathworks.com>...
>
> > Active set with *box constraint* (such as non-negative) requires only matrix-vector product.
> ===============
>
> This is because you can just use the active constraints to eliminate the constrained
> x(i) from the quadratic cost function?
>

No because the projection operator can be trivially computed without any inversion of matrix (min/max is carried out using lower and upper bounds).

Bruno

"Matt " <xys@whatever.com> wrote in message <hdv5s7$csm$1@fred.mathworks.com>...

> W=sum(A.^2,2); %Hessian diagonal
> W=1./W; %invert
>
>
>
> Curvs=A'*(sum(A,2).*W); %majorizing curvatures

This last line should be

Curvs=abs(A)'*(sum(abs(A),2).*W); %majorizing curvatures

"Chad " <parishcm@ornl.gov> wrote in message <hdv1lm$iqp$1@fred.mathworks.com>...
>
> >
> > It might be worth telling us what you want to do with XtX. It's very unusual to need to compute such a matrix explicitly. For example, if all you want is some matrix-vector product y=XtX*b, then you would certainly compute this according to
> >
>
> I often find myself computing X'*X, followed by eigenvectors via eigs(XtX), to do PCA. Is there a way to accomplish PCA without X'*X?

EIGS accepts function handle matrix-vector product, and you don't need to compute X'*X:

X=sprand(100,100,0.1);
Afun=@(b) X'*(X*b)
eigs(Afun,100)

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hdv6s1$e5o$1@fred.mathworks.com>...

> EIGS accepts function handle matrix-vector product, and you don't need to compute X'*X:
>
> X=sprand(100,100,0.1);
> Afun=@(b) X'*(X*b)
> eigs(Afun,100)
>

Strange that SVDS doesn't accept the same. But presumably, since svds is the same as eigs applied to B = [sparse(m,m) A; A' sparse(n,n)], it would be better to
do something like

 X=sprand(100,100,0.1);
 Afun=@(b) [X*b(end+1-size(X,2):end) ; b(1:size(X,1))*X];
 E=eigs(Afun,100);
E=E(E>0);

This way, I'd imagine, you avoid the poorer conditioning of X'*X as opposed to X

"Matt " <xys@whatever.com> wrote in message <hdv8f6$o0q$1@fred.mathworks.com>...

>
> This way, I'd imagine, you avoid the poorer conditioning of X'*X as opposed to X

Matt, many people seem to mix between conditioning number for inversion and eigenvalue problems. They are not necessary related. A matrix has large condition number could be easier to compute for eigenvalues than the opposite. For iterative method, it's sometime preferable to increase condition number, create a larger gap between two neighboring eigenvalues and thus result a faster convergence.

That having said, we note to merit to bring up two different ways of accomplish the calculation of singular values.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hdv9ge$rgp$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <hdv8f6$o0q$1@fred.mathworks.com>...
>
> >
> > This way, I'd imagine, you avoid the poorer conditioning of X'*X as opposed to X
>
> Matt, many people seem to mix between conditioning number for inversion and eigenvalue problems. They are not necessary related. A matrix has large condition number could be easier to compute for eigenvalues than the opposite. For iterative method, it's sometime preferable to increase condition number, create a larger gap between two neighboring eigenvalues and thus result a faster convergence.
==========

But since you can use the eigendecomposition to solve an inversion problem, isn't that a contradiction?

I'll admit that I don't have my head fully wrapped around this, but if I can get a stable eigendecomposition, in spite of a large condition number, then it seems to follow that I can use that decomposition to stably compute an inversion, likewise in spite of the large condition number.

> Note that:
>
> Active set with *box constraint* (such as non-negative) requires only matrix-vector product.
>
> Active set without box constraints (general linear constraints) can be transformed into box-constraint using dual variable. But usually the sparsity of lost because the inverse of the Hessian is needed.
>
> Bruno

so you contradict this reference then:

http://www.cs.utexas.edu/ftp/pub/techreports/tr06-54.pdf

and if i understand you, then they are wrong when they state in section 1.1 (third paragraph) that state-of-the-art fnnls algorithms generally require the computation of A'A? i'm not a "linear algebraist" by any means, but i have implemented just about every piece of freely available source code out there to do nnls / fnnls for my problem, and none come close to the speed/efficiency that Bro's fnnls (and slight variations on it) give me for my nnls problem. every time someone writes a piece of code and says it's better/faster/more efficient, i take it, comment out my call to fnnls, insert their call, then run my program and my computer either locks up or takes at least an order of magnitude longer than Bro's code (which requires A'A).... so maybe you're right and I don't need A'A, but you can see why (from experience) i am skeptical of your claims.... or perhaps i don't need A'A, but
the other methods that save time without calculating A'A are more inefficient in other areas???

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdvk9e$o79$1@fred.mathworks.com>...
> > Note that:
> >
> > Active set with *box constraint* (such as non-negative) requires only matrix-vector product.
> >
> > Active set without box constraints (general linear constraints) can be transformed into box-constraint using dual variable. But usually the sparsity of lost because the inverse of the Hessian is needed.
> >
> > Bruno
>
> so you contradict this reference then:
>
> http://www.cs.utexas.edu/ftp/pub/techreports/tr06-54.pdf
>
> and if i understand you, then they are wrong when they state in section 1.1 (third paragraph) that state-of-the-art fnnls algorithms generally require the computation of A'A?

No I'm not contradict with the above paper, transforming to dual require the need to compute the product of C*inv(A.'*A)*C'*lambda for an given vector of dual lambda, C being the matrix of constraint. An assumption is that this matrix must be strictly positive as well. This is well known (see any good optimization book), but the price is high because the inversion involved. Obviously the article does not consider this transformation, for one of the above inconveniences.

I just try to give you a tip you might use, it's up to use to balance the advantage and inconvenience of such trick in your problem.

Bruno

"Matt " <xys@whatever.com> wrote in message <hdvd2q$704$1@fred.mathworks.com>...

>
> I'll admit that I don't have my head fully wrapped around this, but if I can get a stable eigendecomposition, in spite of a large condition number, then it seems to follow that I can use that decomposition to stably compute an inversion, likewise in spite of the large condition number.

Not exactly, because when you try to bound the norm of the inversion error, the condition number will come to play regardless how well you might know the eigen-decomposition.

Bruno

>
> svds would have been my suggestion, but I'm surprised that resorting to X'*X doesn't give you numerical problems...

Agreed, but it works quite well despite that expectation. My data is non-negative, and can range from somewhat sparse to full (depending on experimental signal/noise). Might these help the stability?

"Chad " <parishcm@ornl.gov> wrote in message <he0qgf$hhl$1@fred.mathworks.com>...
>
> >
> > svds would have been my suggestion, but I'm surprised that resorting to X'*X doesn't give you numerical problems...
>
> Agreed, but it works quite well despite that expectation. My data is non-negative, and can range from somewhat sparse to full (depending on experimental signal/noise). Might these help the stability?
===========

I wouldn't expect non-negativity to help, but sparsity would.

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <hdvk9e$o79$1@fred.mathworks.com>...

every time someone writes a piece of code and says it's better/faster/more efficient, i take it, comment out my call to fnnls, insert their call, then run my program
===================

Did you do this for the algorithm I proposed? I'm curious to know how it compares...

"Matt " <xys@whatever.com> wrote in message <hdv6j6$qst$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <hdv5s7$csm$1@fred.mathworks.com>...
>
> > W=sum(A.^2,2); %Hessian diagonal
> > W=1./W; %invert
> >
> >
> >
> > Curvs=A'*(sum(A,2).*W); %majorizing curvatures
>
> This last line should be
>
> Curvs=abs(A)'*(sum(abs(A),2).*W); %majorizing curvatures

Matt - If you're referring to your code that I reference (in part) above, then yes, I did try that. Sorry to report that my computer grinds to a halt with your code, and I killed the process after about 10 minutes. My 'A' matrices are real, sparse, and on the order of tens or hundreds of thousands of rows by thousands of columns (despite the sparsity, I obviously have lots of RAM and a 64-bit machine!). If I pre-compute A'A and A'b and send it to Bro's fnnls, no problem and it finishes very quickly (well under 1 minute for same problem that I killed after waiting 10 minutes) once that pre-computation is done.

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <he19e4$rfu$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <hdv6j6$qst$1@fred.mathworks.com>...
> > "Matt " <xys@whatever.com> wrote in message <hdv5s7$csm$1@fred.mathworks.com>...
> >
> > > W=sum(A.^2,2); %Hessian diagonal
> > > W=1./W; %invert
> > >
> > >
> > >
> > > Curvs=A'*(sum(A,2).*W); %majorizing curvatures
> >
> > This last line should be
> >
> > Curvs=abs(A)'*(sum(abs(A),2).*W); %majorizing curvatures
>
> Matt - If you're referring to your code that I reference (in part) above, then yes, I did try that. Sorry to report that my computer grinds to a halt with your code, and I killed the process after about 10 minutes.
======

That's very strange. The computations in each iteration of my algorithm are far less intensive than yours, if indeed you are computing A'*A explicitly.
You incorporated the typo fix that I mentioned above?
 

> ======
>
> That's very strange. The computations in each iteration of my algorithm are far less intensive than yours, if indeed you are computing A'*A explicitly.
> You incorporated the typo fix that I mentioned above?
>

Yes, that fix was incorporated. Maybe MATLAB is quite efficient in calculating A'*A if A is sparse enough?

"J RP" <philNOSPAMsonj@hotSPAMMENOTmail.com> wrote in message <he1ddj$6bv$1@fred.mathworks.com>...
> > ======
> >
> > That's very strange. The computations in each iteration of my algorithm are far less intensive than yours, if indeed you are computing A'*A explicitly.
> > You incorporated the typo fix that I mentioned above?
> >
>
> Yes, that fix was incorporated. Maybe MATLAB is quite efficient in calculating A'*A if A is sparse enough?
========================

No, if you have enough computer power to obtain A'*A, then the computations that my proposed algorithm calls for should be trivial.

My impression is that you're not even getting to the main for-loop in my code. Possibly, A and A.' are lying around in your workspace pre-computed, so that when my code calls for you to compute abs(A) (which is also a super large matrix), it starts to strain memory resources.

Below, I've reworked my algorithm somewhat to conserve memory. Whether you continue to play with it is entirely up to you. Just realize that a fair competition between my algorithm and yours would require you to free up enough memory (at minimum by clearing A.') so that the iteration loop is actually reached!

absA=spfun(@(x) abs(x), A); %might be cheaper way to compute abs(A)
Curvatures=( sum(absA,2).' * absA).';
clear absA;

 x=x0; %Initialize x with some vector x0>=0

disp 'Iterations starting now'

for i=1:NumIterations %MAIN LOOP

 gradient=((A*x-b).'*A).';
 x=x-gradient./Curvatures;
 x(x<0)=0;

end

i think you may be correct on the memory issue--i'll keep playing with it, thanks.

> No, if you have enough computer power to obtain A'*A, then the computations that my proposed algorithm calls for should be trivial.
>
> My impression is that you're not even getting to the main for-loop in my code. Possibly, A and A.' are lying around in your workspace pre-computed, so that when my code calls for you to compute abs(A) (which is also a super large matrix), it starts to strain memory resources.
>
> Below, I've reworked my algorithm somewhat to conserve memory. Whether you continue to play with it is entirely up to you. Just realize that a fair competition between my algorithm and yours would require you to free up enough memory (at minimum by clearing A.') so that the iteration loop is actually reached!
>
> absA=spfun(@(x) abs(x), A); %might be cheaper way to compute abs(A)
> Curvatures=( sum(absA,2).' * absA).';
> clear absA;
>
> x=x0; %Initialize x with some vector x0>=0
>
> disp 'Iterations starting now'
>
> for i=1:NumIterations %MAIN LOOP
>
> gradient=((A*x-b).'*A).';
> x=x-gradient./Curvatures;
> x(x<0)=0;
>
> end