Re: bsxfun and sparse matrices

[Jordi Gutiérrez Hermoso]

On 19 May 2012 14:25, Luca Citi <address@hidden> wrote:
> Hello Jordi,
> I was using octave3.2.4 and the following gave me an error:
>
> X = speye(6);
> W = 1:6;
> XW = bsxfun(@times, X, W);
> error: concatenation operator not implemented for ....
>
> Now the HEAD version does not raise the error and works fine.

Btw, bit of terminology: Mercurial doesn't have shouty HEADs. That's a
git and I think svn thing. We have several heads all at once, so I
guess you mean the current head on the default branch. Sometimes, very
rarely, we have more than one head per branch.

In the future, it's best to say which revision you're talking about.
You can see which currently checked-out revision you have with "hg
id", or you can get a bit more verbose output with "hg summary". All
hg commands can be also be abbreviated to their minimal unique initial
substring, so "hg sum" also works.

> The only thing is that octave does not appear to use an heuristic to
> understand if the result can be efficiently saved as sparse or not. In
> fact:
> issparse(XW)
> 0

What does Matlab do here? I'd love to produce sparse output here, but
if Matlab doesn't do it, we're stuck.

> I would like to improve this but if you have time and think this is
> something worth investing time, feel free to do it.

I'm not particularly motivated at present to do this. There's a simple
workaround for this particular case:

    XW = X*diag(W)

I am not sure if this also produces a sparse result in Matlab, since
the sparse result in Octave here depends on our diagonal type, which
doesn't exist in Matlab.

> I am pretty sure it would be more efficient because you know where
> to look in the source code. Otherwise, maybe you can give some
> pointer or suggestions on how to proceed.

Well, the source files you will have to look at are

    src/DLD-FUNCTIONS/bsxfun.cc
    liboctave/bsxfun-defs.cc
    liboctave/bsxfun.h
    liboctave/bsxfun-decl.h

I think the problem is that the DLD bsxfun function doesn't have any
logic for sparse matrices at all. For example, using it on two sparse
matrices also produces a full matrix.

For sparse operators, if you want to make automatic broadcasting work
for them (like how rand(3,3) + rand(3,1) works), we'd have to fix
liboctave/Sparse-op-defs.h first. It's over 2000 lines of pure macro
monstrosity. We don't have to kill its macros before we can work with
it, but I find macros very difficult to work with because I can't run
them in a debugger and I can't get good compiler diagnostics if there
are problems. If you know how to kill macros, I would be very happy to
apply your patches.

HTH,
- Jordi G. H.