![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Hi, there was a thread about this before, diag() is currently only partly useful if you work with numpy-matrices, because the 1d->2d direction doesn't work, as there are no 1d-numpy-matrices. This is unfortunate because a numpy-matrix with shape (n,1) or (1,m) should be naturally treated as a vector, imho. So it would be nice if this could be fixed. It's probably not the most efficient solution, but what I want for numpy-matrix input x is to get: mat(diag(x.A.squeeze)) where diag is the current implementation. This means that if x is not a vector ("truly 2d"), then nothing is changed. But if one of the dimensions of x is ==1, then it's turned into a 1d-array, and diag works as it should. Does that sound reasonable? Thanks, Sven
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Sven Schreiber wrote:
Not for numpy.diag() in my opinion. However, I won't object to a numpy.matlib.diag() that knows about matrix objects and behaves the way you want. -- 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
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Robert Kern schrieb:
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 repeat, the matlib solution would be ok for me, but in some sense not fixing numpy.diag could contribute to the feeling of matrices being only second-class citizens. cheers, Sven
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Sven Schreiber wrote:
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. 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. -- 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
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Robert Kern schrieb:
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 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... cheers, Sven
![](https://secure.gravatar.com/avatar/6a1dc50b8d79fe3b9a5e9f5d8a118901.jpg?s=120&d=mm&r=g)
Every time I want to make a diagonal matrix from a vector I have to do a google search. So I'm with you Sven: diag(NxN matrix) should return a Nx1 matrix diag(Nx1 or 1xN matrix) should return a NxN matrix instead of the current behavior: diag(NxN matrix) returns an array diag(Nx1 matrix) returns a 1x1 array
x
matrix([[ 0.43298158, 0.84572719], [ 0.1391546 , 0.66412104]])
On 7/28/06, Sven Schreiber <svetosch@gmx.net> wrote:
![](https://secure.gravatar.com/avatar/5dde29b54a3f1b76b2541d0a4a9b232c.jpg?s=120&d=mm&r=g)
Keith Goodman wrote:
diag(NxN matrix) should return a Nx1 matrix diag(Nx1 or 1xN matrix) should return a NxN matrix
This is the key problem: extracting the diagonal of a matrix and creating a matrix from a diagonal are two different operations: overloading one function to do both was a bad idea to begin with. Maybe we should just keep diag() as is is for backward compatibility (deprecated), and make: get_diag() and make_diag() instead. Then it would be unambiguous what you wanted with: make_diag(<Nx1array>) You can call them something else, but you get the idea. -Chris -- Christopher Barker, Ph.D. Oceanographer NOAA/OR&R/HAZMAT (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Christopher Barker wrote:
As .diagonal() is a method on both arrays and matrices, we don't need a separate get_diag(). In the abstract, I heartily approve. However, I'm not sure that deprecation in favor of a clumsier name is going to be effective in practice. diag([a, b, c]) is just too tempting and will remain an attractive nuisance. It matches mathematical notation ever so nicely. -- 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
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Here's my attempt at summarizing the diag-discussion. The symptom is that currently transforming something like the vector a b c into the diagonal matrix a 0 0 0 b 0 0 0 c is not directly possible if you're working with numpy-matrices (i.e. the vector is of type matrix and has shape (n,1) or (1,n)). Blindly doing a diag(vector) just gives you the first element of the vector. Instead you need to do something like mat(diag(vector.A.squeeze). Proposals that were put forward: 1) Modify the existing diag() to treat not only 1d-arrays as vectors, but also numpy-matrices of shapes (1,n) and (n,1), and therefore have diag() turn them into a diagonal matrix with all their elements. This would be my preferred solution. Robert is against this, because it is not 100% backwards compatible and because it requires the diag-function to differentiate between the input types. 2) Deprecate the use of diag which is overloaded with making diagonal matrices as well as getting diagonals. Instead, use the existing .diagonal() for getting a diagonal, and introduce a new make_diag() function which could easily work for numpy-arrays and numpy-matrices alike. This is a blend of Christopher's and Robert's comments. As robert noted, the problem is that diag() will probably be more attractive than make_diag() for numpy-array users, and so make_diag() would effectively only be used by numpy-matrix users, a split which would not be nice. An effectively equivalent solution would then be a new method .make_diag() for numpy-matrices, which was suggested by Robert (if I understood correctly). 3) Add a new function diag() to numpy.matlib implementing the workaround, something like: def diag(vector): return mat(N.diag(vector.A.squeeze) (not sure if all namespace issues are avoided like that, or if it's efficient) That way, numpy.matlib.diag() could be used with numpy-matrices just like numpy.diag() with numpy-arrays. Side effect: this even turns a 1d-array into a numpy-matrix, which is useful e.g. for dealing with eigenvalues (which are 1d-arrays even for numpy-matrix input). The third solution is probably the way of least resistance. However, afaik this would be the first time that numpy needs different functions to perform essentially the same operation. Since Robert already cited the obligatory zen phrase, what about this special casing? Anyway, any of those solutions is better than the status quo imo. I'd be happy if one of them could be implemented. When I get back from my vacation and in case nothing has been changed, I'm going to file a ticket so that it doesn't get forgotten. Thanks for your patience and have a good time! Sven Robert Kern schrieb:
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
Hi Sven, On 7/28/06, Sven Schreiber <svetosch@gmx.net> wrote:
Here's my attempt at summarizing the diag-discussion.
<snip> 2) Deprecate the use of diag which is overloaded with making diagonal
This would be my preference, but with functions {get,put}diag. We could also add a method or function asdiag, which would always return a diagonal matrix made from *all* the elements of the matrix taken in order. For (1,n) or (n,1) this would do what you want. For other matrices the result would be something new and probably useless, but at least it wouldn't hurt. Chuck
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Charles R Harris schrieb:
This seems to have been implemented now by the new diagflat() function. So, matrix users can now use m.diagonal() for the matrix->vector direction of diag(), and diagflat(v) for the vector->matrix side of diag(), and always get numpy-matrix output for numpy-matrix input. Thanks a lot for making this possible! One (really minor) comment: "diagflat" as a name is not optimal imho. Are other suggestions welcome, or is there a compelling reason for this name? Thanks, sven
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
The following are from pyGAUSS. fwiw, Alan Isaac #diag: diagonal of matrix as column vector (2D only!) def diag(x): if not isinstance(x,numpy.matrix): x = numpy.asanyarray(x) assert(len(x.shape)==2), "For 2-d arrays only." return x.diagonal(offset=0,axis1=-2,axis2=-1).reshape((-1,1)) else: return x.diagonal().T #diagrv: insert v as diagonal of matrix x (2D only!) def diagrv(x,v,copy=True): x = numpy.matrix( x, copy=copy ) stride = 1 + x.shape[1] x.flat[ slice(0,x.size,stride) ] = v return x
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Sven Schreiber wrote:
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.
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
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Sven Schreiber wrote:
Not for numpy.diag() in my opinion. However, I won't object to a numpy.matlib.diag() that knows about matrix objects and behaves the way you want. -- 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
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Robert Kern schrieb:
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 repeat, the matlib solution would be ok for me, but in some sense not fixing numpy.diag could contribute to the feeling of matrices being only second-class citizens. cheers, Sven
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Sven Schreiber wrote:
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. 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. -- 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
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Robert Kern schrieb:
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 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... cheers, Sven
![](https://secure.gravatar.com/avatar/6a1dc50b8d79fe3b9a5e9f5d8a118901.jpg?s=120&d=mm&r=g)
Every time I want to make a diagonal matrix from a vector I have to do a google search. So I'm with you Sven: diag(NxN matrix) should return a Nx1 matrix diag(Nx1 or 1xN matrix) should return a NxN matrix instead of the current behavior: diag(NxN matrix) returns an array diag(Nx1 matrix) returns a 1x1 array
x
matrix([[ 0.43298158, 0.84572719], [ 0.1391546 , 0.66412104]])
On 7/28/06, Sven Schreiber <svetosch@gmx.net> wrote:
![](https://secure.gravatar.com/avatar/5dde29b54a3f1b76b2541d0a4a9b232c.jpg?s=120&d=mm&r=g)
Keith Goodman wrote:
diag(NxN matrix) should return a Nx1 matrix diag(Nx1 or 1xN matrix) should return a NxN matrix
This is the key problem: extracting the diagonal of a matrix and creating a matrix from a diagonal are two different operations: overloading one function to do both was a bad idea to begin with. Maybe we should just keep diag() as is is for backward compatibility (deprecated), and make: get_diag() and make_diag() instead. Then it would be unambiguous what you wanted with: make_diag(<Nx1array>) You can call them something else, but you get the idea. -Chris -- Christopher Barker, Ph.D. Oceanographer NOAA/OR&R/HAZMAT (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Christopher Barker wrote:
As .diagonal() is a method on both arrays and matrices, we don't need a separate get_diag(). In the abstract, I heartily approve. However, I'm not sure that deprecation in favor of a clumsier name is going to be effective in practice. diag([a, b, c]) is just too tempting and will remain an attractive nuisance. It matches mathematical notation ever so nicely. -- 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
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Here's my attempt at summarizing the diag-discussion. The symptom is that currently transforming something like the vector a b c into the diagonal matrix a 0 0 0 b 0 0 0 c is not directly possible if you're working with numpy-matrices (i.e. the vector is of type matrix and has shape (n,1) or (1,n)). Blindly doing a diag(vector) just gives you the first element of the vector. Instead you need to do something like mat(diag(vector.A.squeeze). Proposals that were put forward: 1) Modify the existing diag() to treat not only 1d-arrays as vectors, but also numpy-matrices of shapes (1,n) and (n,1), and therefore have diag() turn them into a diagonal matrix with all their elements. This would be my preferred solution. Robert is against this, because it is not 100% backwards compatible and because it requires the diag-function to differentiate between the input types. 2) Deprecate the use of diag which is overloaded with making diagonal matrices as well as getting diagonals. Instead, use the existing .diagonal() for getting a diagonal, and introduce a new make_diag() function which could easily work for numpy-arrays and numpy-matrices alike. This is a blend of Christopher's and Robert's comments. As robert noted, the problem is that diag() will probably be more attractive than make_diag() for numpy-array users, and so make_diag() would effectively only be used by numpy-matrix users, a split which would not be nice. An effectively equivalent solution would then be a new method .make_diag() for numpy-matrices, which was suggested by Robert (if I understood correctly). 3) Add a new function diag() to numpy.matlib implementing the workaround, something like: def diag(vector): return mat(N.diag(vector.A.squeeze) (not sure if all namespace issues are avoided like that, or if it's efficient) That way, numpy.matlib.diag() could be used with numpy-matrices just like numpy.diag() with numpy-arrays. Side effect: this even turns a 1d-array into a numpy-matrix, which is useful e.g. for dealing with eigenvalues (which are 1d-arrays even for numpy-matrix input). The third solution is probably the way of least resistance. However, afaik this would be the first time that numpy needs different functions to perform essentially the same operation. Since Robert already cited the obligatory zen phrase, what about this special casing? Anyway, any of those solutions is better than the status quo imo. I'd be happy if one of them could be implemented. When I get back from my vacation and in case nothing has been changed, I'm going to file a ticket so that it doesn't get forgotten. Thanks for your patience and have a good time! Sven Robert Kern schrieb:
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
Hi Sven, On 7/28/06, Sven Schreiber <svetosch@gmx.net> wrote:
Here's my attempt at summarizing the diag-discussion.
<snip> 2) Deprecate the use of diag which is overloaded with making diagonal
This would be my preference, but with functions {get,put}diag. We could also add a method or function asdiag, which would always return a diagonal matrix made from *all* the elements of the matrix taken in order. For (1,n) or (n,1) this would do what you want. For other matrices the result would be something new and probably useless, but at least it wouldn't hurt. Chuck
![](https://secure.gravatar.com/avatar/2b8a71c6c4c0df90de065ce45ee9df33.jpg?s=120&d=mm&r=g)
Charles R Harris schrieb:
This seems to have been implemented now by the new diagflat() function. So, matrix users can now use m.diagonal() for the matrix->vector direction of diag(), and diagflat(v) for the vector->matrix side of diag(), and always get numpy-matrix output for numpy-matrix input. Thanks a lot for making this possible! One (really minor) comment: "diagflat" as a name is not optimal imho. Are other suggestions welcome, or is there a compelling reason for this name? Thanks, sven
![](https://secure.gravatar.com/avatar/39916bae984cb93b797efd2b175f59c0.jpg?s=120&d=mm&r=g)
The following are from pyGAUSS. fwiw, Alan Isaac #diag: diagonal of matrix as column vector (2D only!) def diag(x): if not isinstance(x,numpy.matrix): x = numpy.asanyarray(x) assert(len(x.shape)==2), "For 2-d arrays only." return x.diagonal(offset=0,axis1=-2,axis2=-1).reshape((-1,1)) else: return x.diagonal().T #diagrv: insert v as diagonal of matrix x (2D only!) def diagrv(x,v,copy=True): x = numpy.matrix( x, copy=copy ) stride = 1 + x.shape[1] x.flat[ slice(0,x.size,stride) ] = v return x
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Sven Schreiber wrote:
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.
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
participants (6)
-
Alan G Isaac
-
Charles R Harris
-
Christopher Barker
-
Keith Goodman
-
Robert Kern
-
Sven Schreiber