[Numpy-discussion] fixing diag() for matrices

Sven Schreiber svetosch at gmx.net
Fri Jul 28 07:46:14 CDT 2006

Robert Kern schrieb:
> Sven Schreiber wrote:
>> That would be fine with me. However, I'd like to point out that after
>> some bug-squashing currently all numpy functions deal with
>> numpy-matrices correctly, afaik. The current behavior of numpy.diag
>> could be viewed as a violation of that principle. (Because if x has
>> shape (n,1), diag(x) returns only the first entry, which is pretty
>> stupid for a diag-function operating on a vector.)
> I don't think so. It's operating exactly as it is documented to. It doesn't do 
> what you want, but it's not supposed to read your mind and guess that you are 
> treating (n,1) and (1,n) arrays as linear algebraic vectors and that you want to 
> form a diagonal matrix from that vector. It handles matrix objects just fine; it 
> does not obey a particular convention that you want to use, though. That's 
> neither broken nor a violation of the principle that most numpy functions should 
> accept matrix objects; it's just not what you want.

Ok, let me get some facts straight:

1. If you're using numpy-matrices, algebraic vectors are *automatically*
(n,1) or (1,n). I didn't make it up, and nobody has to read my mind to
understand that; it's just a *fact* of numpy (and you knew that long
before I was even aware of numpy's existence ;-).

2. I never claimed that the documentation of diag is wrong. I'm just
questioning whether its behavior with vectors represented as
numpy-matrices (which have shape (n,1) or (1,n) and are therefore always
2d in the numpy sense) is really intended, because I currently doubt
that it's useful for *anyone*. (Of course you can prove me wrong there...)

3. The cited principle is not (only) about accepting matrices, it's
about "honoring" them by either preserving their type for the output,
and of course by doing the equivalent thing as for the equivalent
numpy-array input. Currently, however, diag() is not providing the
"vector-to-diagonal-matrix" functionality if you work with
numpy-matrices (modulo some obvious workarounds). To me, that's a partly
broken function, and it's *not* handling matrix objects "just fine".

> I don't want to introduce a backwards-compatibility-breaking special case to the 
> function. "Special cases aren't special enough to break the rules." Different 
> functionality should go into a different function.

I hope (but somehow doubt...) that I've convinced you that it's actually
about applying the same logical rule to all input types, not about
introducing a special case. I agree that in principle backwards
compatibility could be an issue; however that would mean that people are
actually using the strange behavior that diag() displays with (n,1) or
(1,n) matrix-vectors (namely returning only the first element).

Does anybody do that? I doubt it, but of course I can't prove it; does
anybody on the list know of an example where it's used?

Having said all that in this lengthy message, I don't want to push it
too far. I believe that the right thing to do would be to fix diag()
along the lines I suggested. If I haven't managed to convince you and/or
the other developers, so be it, and I can live with a new
numpy.matlib.diag() just fine.

In the latter case, may I also suggest that a new numpy.matlib.diag()
always return a numpy-matrix even when given a numpy-array? Background:
Currently (and I think it's ok that way) the eig*() functions return the
eigenvalues in a 1d-array, even for numpy-matrices as inputs. It is
fairly usual to use a diagonal matrix with the eigenvalues on the
diagonal in an algebraic formula. That could be achieved with
n.matlib.diag(n.linalg.eig*(...)[0]) without an enclosing mat() call
only if n.matlib.diag always returns a numpy-matrix.

Thanks for reading until the end...

More information about the Numpy-discussion mailing list