[Numpy-discussion] fixing diag() for matrices

[Robert Kern]

Sven Schreiber wrote:
> 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".

Unfortunately, there is a many-to-one correspondence between numpy-array inputs 
and matrix-array inputs. Going the other way is ambiguous. There is still some 
guessing involved.

>> 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.

But you can't code it without special casing. None of the functions should be 
written like this:

   if type(input) is matrix:
     # do this one thing
   else:
     # do a different thing

They belong in two different functions. Or methods. You can't extend it to the n 
different subclasses of ndarray that people might want to use. This kind of 
special casing ought to give you screaming jibblies. The current behavior of 
diag() already gives me the screaming jibblies given that it both extracts a 
diagonal from an MxN array and creates a NxN array from an N array. <shudder>

If I wrote an API like that at work, my project manager would come over to my 
office and Have Words with me. And then tell me to rewrite it. It's simply bad 
design.

I think making a method on matrix objects that does what you want is probably 
your best bet.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco