Path: news.mathworks.com!not-for-mail
From: "Tim Davis" <davis@cise.ufl.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Sparse matrix product again
Date: Fri, 2 May 2008 15:47:03 +0000 (UTC)
Organization: University of Florida
Lines: 82
Message-ID: <fvfd1n$1dc$1@fred.mathworks.com>
References: <67rg20F2pn19mU1@mid.dfncis.de>
Reply-To: "Tim Davis" <davis@cise.ufl.edu>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1209743223 1452 172.30.248.35 (2 May 2008 15:47:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 2 May 2008 15:47:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 45902
Xref: news.mathworks.com comp.soft-sys.matlab:466295


Joachim Selke <j.selke@tu-bs.de> wrote in message
<67rg20F2pn19mU1@mid.dfncis.de>...
> Hi,
> 
> in <news:frk8n2$dqp$1@newsserver.rrzn.uni-hannover.de> I
asked how to
> efficiently compute a "sparse matrix product":
> 
> Given three matrices
> 
>   A: m x k   (full)
>   B: k x n   (full)
>   X: m x n   (sparse, around 2% density)
> 
> where m and n are very large numbers and k is small, I
want to compute
> the matrix product of A and B but I only need the entries
of A * B at
> the places where X is non-zero.
> 
> Thanks to Tim Davis and Roger Stafford I got two nice
solutions that I
> combined and optimized a bit further. My current best
solution is the
> following ("best" means fastest and memory-efficient enough):
> 
> % Basic idea: Construct the sparse product matrix row by row.
> 
> Atr = A';
> [rows, cols, vals] = find(X);
> lX = logical(X);
> nnz_bycol = full(sum(lX));
> cumnnz_bycol = [0 cumsum(nnz_bycol)];
> 
> for j = 1:n
>     t_start = cumnnz_bycol(j) + 1;
>     t_end = cumnnz_bycol(j + 1);
>     rows_j = rows(t_start:t_end);
>     vals(t_start:t_end) = Atr(:, rows_j)' * B(:, j);
> end
> 
> AB = sparse(rows, cols, vals, m, n);
> 
> 
> Now I did some experiments using random matrices (m =
50000, n = 10000,
> k = 15). The total time needed for constructing the
product was around
> 2.3 seconds on my computer. When looking at the details, I
found that
> almost the complete computation time is used by the statement
> "Atr(:, rows_j)." So, obviously this algorithm is extremely
> memory-bound. I tried a lot but was not able to solve this
problem.
> 
> Are there any known tricks what can I do to speed-up the
memory accesses?
> 
> (I should note that I want to compute this product for
many different
> matrices A and B while keeping X fixed. Maybe this fact
can be used in
> optimization ...)
> 
> 
> Joachim

That looks about as good as you could get, to me.

The Atr(:,rows_j)'*B(:j) is a neat trick; using a
matrix-vector multiply to do lots of entries in AB at once.
 I can't think of a better way ... except to go to a
mexFunction of course.

The reason I say that is because a mexFunction could do
Atr(:,rows_j)'*B(:j) without forming the Atr(:,rows_j)
matrix as an intermediate step.

The subsref and subsasgn operations in MATLAB aren't speed
demons, unfortunately, particularly for the sparse case
(this is the dense case, though, since Atr is dense).