Thread: [Numpy-discussion] How do I make a diagonal matrix?

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 422-6466

How do I make a NxN diagonal matrix with a Nx1 column vector x along
the diagonal?

Keith Goodman schrieb:
> How do I make a NxN diagonal matrix with a Nx1 column vector x along
> the diagonal?
> 

>>> help(n.diag)
Help on function diag in module numpy.lib.twodim_base:

diag(v, k=0)
    returns the k-th diagonal if v is a array or returns a array
    with v as the k-th diagonal if v is a vector.

On Friday 23 June 2006 16:34, Sven Schreiber wrote:
   [SS]: Keith Goodman schrieb:
   [SS]: > How do I make a NxN diagonal matrix with a Nx1 column vector x along
   [SS]: > the diagonal?
   [SS]: > 
   [SS]: 
   [SS]: >>> help(n.diag)
   [SS]: Help on function diag in module numpy.lib.twodim_base:
   [SS]: 
   [SS]: diag(v, k=0)
   [SS]:     returns the k-th diagonal if v is a array or returns a array
   [SS]:     with v as the k-th diagonal if v is a vector.

See also the Numpy Example List for a few examples:
http://www.scipy.org/Numpy_Example_List

J.

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

On 6/23/06, Sven Schreiber <sve...@gm...> wrote:
> Keith Goodman schrieb:
> > How do I make a NxN diagonal matrix with a Nx1 column vector x along
> > the diagonal?
> >
>
> >>> help(n.diag)
> Help on function diag in module numpy.lib.twodim_base:
>
> diag(v, k=0)
>     returns the k-th diagonal if v is a array or returns a array
>     with v as the k-th diagonal if v is a vector.

I tried

>> x = rand(3,1)

>> diag(x)
array([ 0.87113114])

Isn't rand(3,1) a vector?

Off list I was given the example:
 x=rand(3)
diag(3)

That works. But my x is a Nx1 matrix. I can't get it to work with matrices.

Joris: The Numpy Example List looks good. I hadn't come across that before.

Keith Goodman schrieb:

> 
> Isn't rand(3,1) a vector?

afaik not in numpy's terms, because two numbers are given for the
dimensions -- I also struggle with that, because I'm a matrix guy like
you ;-)

> 
> Off list I was given the example:
>  x=rand(3)
> diag(3)
> 
> That works. But my x is a Nx1 matrix. I can't get it to work with matrices.
> 

ok, good point; with your x then diag(x.A[:,0]) should work, although
it's not very pretty. Maybe there are better ways, but I agree it would
be nice to be able to use matrices directly.

-sven

On Fri, 23 Jun 2006, Keith Goodman apparently wrote:=20
> my x is a Nx1 matrix. I can't get it to work with matrices.=20

Hmm. One would think that diag() would accept a flatiter=20
object, but it does not.  Shouldn't it??

But anyway, you can squeeze x:

>>> x
matrix([[ 0.46474951],
       [ 0.0688041 ],
       [ 0.61141623]])
>>> y=3DN.diag(N.squeeze(x.A))
>>> y
array([[ 0.46474951,  0.        ,  0.        ],
       [ 0.        ,  0.0688041 ,  0.        ],
       [ 0.        ,  0.        ,  0.61141623]])

hth,
Alan Isaac

On Fri, 23 Jun 2006, Alan G Isaac apparently wrote:=20
> you can squeeze x=20

True, but a silly solution.
Alan

Alan G Isaac wrote:

>On Fri, 23 Jun 2006, Keith Goodman apparently wrote: 
>  
>
>>my x is a Nx1 matrix. I can't get it to work with matrices. 
>>    
>>
>
>Hmm. One would think that diag() would accept a flatiter 
>object, but it does not.  Shouldn't it??
>  
>

It doesn't?  

try:

a = rand(3,4)
diag(a.flat).shape

which prints
(12,12)
for me.

Also:

 >>> a = ones((2,3))

 >>> diag(a.flat)
array([[1, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0],
       [0, 0, 0, 1, 0, 0],
       [0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 1]])

> Alan G Isaac wrote:=20
>> Hmm. One would think that diag() would accept a flatiter=20
>> object, but it does not.  Shouldn't it??=20

On Fri, 23 Jun 2006, Travis Oliphant apparently wrote:=20
> It doesn't?=20
> try:=20
> a =3D rand(3,4)=20
> diag(a.flat).shape=20

OK, but then try:
>>> a=3DN.mat(a)
>>> N.diag(a.flat).shape
(1,)

Why is a.flat not the same as a.A.flat?

Alan Isaac

Alan G Isaac wrote:

>>Alan G Isaac wrote: 
>>    
>>
>>>Hmm. One would think that diag() would accept a flatiter 
>>>object, but it does not.  Shouldn't it?? 
>>>      
>>>
>
>
>On Fri, 23 Jun 2006, Travis Oliphant apparently wrote: 
>  
>
>>It doesn't? 
>>try: 
>>a = rand(3,4) 
>>diag(a.flat).shape 
>>    
>>
>
>OK, but then try:
>  
>
>>>>a=N.mat(a)
>>>>N.diag(a.flat).shape
>>>>        
>>>>
>(1,)
>
>Why is a.flat not the same as a.A.flat?
>  
>

It is the same object except for the pointer to the underlying array.  
When asarray(a.flat) get's called it looks to the underlying array to 
get the sub-class and constructs that sub-class (and matrices can never 
be 1-d).  Thus, it's a "feature"

-Travis

> Alan G Isaac wrote:=20
>> Why is a.flat not the same as a.A.flat?=20

On Fri, 23 Jun 2006, Travis Oliphant apparently wrote:=20
> It is the same object except for the pointer to the=20
> underlying array.  When asarray(a.flat) get's called it=20
> looks to the underlying array to get the sub-class and=20
> constructs that sub-class (and matrices can never be 1-d). =20
> Thus, it's a "feature"=20

I doubt I will prove the only one to stumble over this.

I can roughly understand why a.ravel() returns a matrix;
but is there a good reason to forbid truly flattening the matrix?

My instincts are that a flatiter object should not have this=20
hidden "feature": flatiter objects should produce=20
a consistent behavior in all settings, regardless of the=20
underlying array. Anything else will prove too surprising.

fwiw,
Alan

Alan G Isaac wrote:

>>Alan G Isaac wrote: 
>>    
>>
>>>Why is a.flat not the same as a.A.flat? 
>>>      
>>>
>
>
>On Fri, 23 Jun 2006, Travis Oliphant apparently wrote: 
>  
>
>>It is the same object except for the pointer to the 
>>underlying array.  When asarray(a.flat) get's called it 
>>looks to the underlying array to get the sub-class and 
>>constructs that sub-class (and matrices can never be 1-d).  
>>Thus, it's a "feature" 
>>    
>>
>
>
>I doubt I will prove the only one to stumble over this.
>
>I can roughly understand why a.ravel() returns a matrix;
>but is there a good reason to forbid truly flattening the matrix?
>  
>
Because matrices are never 1-d.  This is actually pretty consistent 
behavior.

>My instincts are that a flatiter object should not have this 
>hidden "feature": flatiter objects should produce 
>a consistent behavior in all settings, regardless of the 
>underlying array. Anything else will prove too surprising.
>  
>
I think you are right that this is a bug, though.  Because __array__() 
(which is where the behavior comes from) should return a base-class 
array (not a sub-class).

-Travis

Travis Oliphant wrote:

>Alan G Isaac wrote:
>  
>
>>    
>>
>I think you are right that this is a bug, though.  Because __array__() 
>(which is where the behavior comes from) should return a base-class 
>array (not a sub-class).
>  
>
This is fixed in SVN.

-Travis

> Alan G Isaac wrote:=20
>> I can roughly understand why a.ravel() returns a matrix;=20
>> but is there a good reason to forbid truly flattening the matrix?=20

On Fri, 23 Jun 2006, Travis Oliphant apparently wrote:=20
> Because matrices are never 1-d.  This is actually pretty=20
> consistent behavior.=20

Yes; that's why I can understand ravel.
But I was referring to flat with the question.

On Fri, 23 Jun 2006, Travis Oliphant apparently wrote:=20
> I think you are right that this is a bug, though.  Because=20
> __array__() (which is where the behavior comes from)=20
> should return a base-class array (not a sub-class).=20

Thanks for fixing this!!

Alan

On Fri, Jun 23, 2006 at 07:55:47AM -0700, Keith Goodman wrote:
> On 6/23/06, Sven Schreiber <sve...@gm...> wrote:
> > Keith Goodman schrieb:
> > > How do I make a NxN diagonal matrix with a Nx1 column vector x along
> > > the diagonal?
> > >
> >
> > >>> help(n.diag)
> > Help on function diag in module numpy.lib.twodim_base:
> >
> > diag(v, k=3D0)
> >     returns the k-th diagonal if v is a array or returns a array
> >     with v as the k-th diagonal if v is a vector.
>=20
> I tried
>=20
> >> x =3D rand(3,1)
>=20
> >> diag(x)
> array([ 0.87113114])
>=20
> Isn't rand(3,1) a vector?

No:
In [13]: rand(3).shape =20
Out[13]: (3,)

In [14]: rand(3,1).shape
Out[14]: (3, 1)

A "vector" is an array with only one dimension. Here, you have a 3x1
"matrix"...

>=20
> Off list I was given the example:
>  x=3Drand(3)
> diag(3)

So you've got the solution!

> That works. But my x is a Nx1 matrix. I can't get it to work with matrice=
s.

???
Don't understand what you cannot make work, here.

In [15]: x=3Drand(3,1)    =20

In [18]: diag(x[:,0])
Out[18]:=20
array([[ 0.2287158 ,  0.        ,  0.        ],
       [ 0.        ,  0.50571537,  0.        ],
       [ 0.        ,  0.        ,  0.72304857]])

What else would you like?

David

> Joris: The Numpy Example List looks good. I hadn't come across that befor=
e.
>=20

David Douard                             LOGILAB, Paris (France)
Formations Python, Zope, Plone, Debian : http://www.logilab.fr/formations
D=E9veloppement logiciel sur mesure :      http://www.logilab.fr/services
Informatique scientifique :              http://www.logilab.fr/science

On Fri, 23 Jun 2006, Sven Schreiber apparently wrote:=20
>>>> help(n.diag)=20
> Help on function diag in module numpy.lib.twodim_base:=20
> diag(v, k=3D0)=20
>     returns the k-th diagonal if v is a array or returns a array=20
>     with v as the k-th diagonal if v is a vector.=20

That is pretty damn obscure.
Apparently Travis's new doc string did not survive?
The Numpy book says:
diag (v, k=3D0)
    Return the kth diagonal if v is a 2-d array, or returns=20
    an array with v as the kth diagonal if v is a 1-d array.

That is better but not great.
I think what is wanted is:
diag (v, k=3D0)
    If v is a 2-d array:
    return a copy of the kth diagonal of v (as a 1-d array).
    If v is a 1-d array:
    return a 2-d array with a copy of v as the kth diagonal=20
    (and zeros elsewhere).

fwiw,
Alan Isaac

PS As a response to the question, it might be worth noting=20
the following.
>>> y=3DN.zeros((5,5))
>>> values=3DN.arange(1,6)
>>> indices=3Dslice(0,25,6)
>>> y.flat[indices]=3Dvalues
>>> y
array([[1, 0, 0, 0, 0],
       [0, 2, 0, 0, 0],
       [0, 0, 3, 0, 0],
       [0, 0, 0, 4, 0],
       [0, 0, 0, 0, 5]])

Generalizing we end up with the following (from pyGAUSS):

def diagrv(x,v,copy=3DTrue):
=09if copy: x =3D numpy.matrix( x, copy=3DTrue )
=09else: x =3D numpy.matrix( x, copy=3DFalse )
=09stride =3D 1 + x.shape[1]
=09x.flat[ slice(0,x.size,stride) ] =3D v
        return x

Thread: [Numpy-discussion] How do I make a diagonal matrix?

A package for scientific computing with Python

numpy-discussion