Re: deprecated functions

[Jaroslav Hajek]

On Wed, Mar 4, 2009 at 11:26 PM, John W. Eaton <address@hidden> wrote:
> On  4-Mar-2009, Jason Riedy wrote:
>
> | And John W. Eaton writes:
> | > Maybe I'm slow, but what is the large impact of having diag x sparse
> | > operations convert to full?  What does it break?  I guess I could use
> | > some examples.
> |
> | A typical use (for me) of diag * sparse is to rescale the matrix before
> | solving a system with it.  Many sparse matrices from automatic systems
> | (e.g. optimization problem builders) are horribly ill-scaled and appear
> | to be ill-conditioned.  Rescaling the problem can greatly improve the
> | solution.  The rescaling should have negligible cost.
> |
> | For a direct method, if diag * sparse is full, you're lead to a full
> | factorization and needlessly slow performance.  For an iterative method,
> | making A dense makes A*x far, far more expensive than it should be.
> |
> | Yes, code *could* contain an if sparse(...) test, but why?  What is the
> | benefit of forcing code to care about structure when it's otherwise
> | agnostic?
>
> I guess my question wasn't clear.  I understand why it would be useful
> to have sparse x diag operations preserve sparsity.  What I'm asking
> is what would be the big problem with doing that?  How does it cause
> significant compatibility problems.  I realize that if we made the
> change so that
>
>  diag (v) * sparse_matrix  -->  sparse in Octave, full in Matlab
>
> that there would be a difference.  But exactly what would be
> catastrophically bad about it?

Well the (possibly) bad thing is that you may get sparse matrices
where you used to get dense ones.
And, consequently, the well known "angry complaints". I don't think
there's a negative impact from any
other perspective, rather a positive one.

> After all, we currently have diag (v) producing a specialized matrix
> type instead of a full matrix, so why is that not also as bad for
> compatibility?  What makes the sparse x diag --> sparse worse?
>

It is not so bad because diagonal matrices will convert themselves
implicitly to full matrices in most occasions, unless the operation is
specialized. They will yield mostly identical results to full matrices
in their place, save from the recently discussed NaN issues (in which
case they return a better result).
For instance, should diag + sparse return a full (compatible) or
sparse (reasonable) result?

That being said, I don't oppose such a change. It would be very
logical, now that we have those special matrices. I'm just warning
that it breaks compatibility. But maybe others will oppose.

The following applies if we agree on this change:
I noticed Jason has already started working on that part. Jason, can
you share the details? We would need mixed diagonal-sparse (*,\,+,-)
and permutation-sparse operations implemented in liboctave, added
probably to Sparse-op-defs.h and sparse-mx-ops.
The diagonal divisors should probably go to src/sparse-xdiv.cc. You
can inspect the divisors in xdiv.cc to see what behavior is needed -
the diagonal matrices, when used for division, should return a
least-squares minimum-norm solution, aka pseudoinverse. That means,
the rectangular part and possible zero diagonal elements need to be
properly ignored, or they can generate zero blocks in the result,
depending on orientation.

As for the interpreter part, I think you should be able to reuse
OPERATORS/op-dm-template.cc and I would strongly advise you to do so.

regards


-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz