[SciPy-user] Easy way to make a block diagonal matrix?

[josef.pktd at gmail.com]

On Thu, May 21, 2009 at 10:18 AM, Ivo Maljevic <ivo.maljevic at gmail.com> wrote:
>> I disagree because block diagonal does have a special meaning and the
>> result is not a diagonal matrix!
>>
>> > If all individual component matrices are square, then you get the
>> > wikipedia definition.
>> >
>> My understanding is that only if the inputs are diagonal matrices will
>> you get a block diagonal matrix from this function.
>>
>
> I am not a scipy developer, but from time to time I send an email to this
> list with either questions or answers to somebody
> else's questions.
>
> What do you mean when you say "My understanding is that only if the inputs
> are diagonal matrices will
> you get a block diagonal matrix from this function."
>
> Here is a non-definition block diagonal matrix with functionality similar to
> matlab:
>
>>>> a=numpy.ones([2,2])
>>>> b=numpy.random.rand(3,3)
>>>> import block
>>>> c=block.block_diag(a,b)
>>>> print c
> [[ 1.          1.          0.          0.          0.        ]
>  [ 1.          1.          0.          0.          0.        ]
>  [ 0.          0.          0.93924665  0.43404552  0.46698808]
>  [ 0.          0.          0.61331601  0.23593332  0.39016641]
>  [ 0.          0.          0.10644194  0.66638397  0.6305998 ]]
>
>
> as you can see, 'a' and 'b' are not diagonal matrices. I think the only
> question is whether this function should work only in square matrix cases
> such as next example:
>
>>>> bb=numpy.random.randn(2,2)
>>>> cc=block.block_diag(a,bb)
>>>> print cc
> [[ 1.          1.          0.          0.        ]
>  [ 1.          1.          0.          0.        ]
>  [ 0.          0.          0.57135707 -0.03777517]
>  [ 0.          0.         -0.22069926 -1.40840111]]
>>>>
>

If the user uses only square matrices, then (s)he will get only "block
diagonal" (according to wikipedia) matrices back.

I don't see why we should choose another name just because the
function allows for more flexibility, and I don't see a use case for
automatically checking whether the matrices are square, so that the
user does not violate the square definition.

Josef

another example

>>> x = block_diag(np.eye(2), [[1,2],[3,4],[5,6]])
>>> x.shape
(5, 4)
>>> np.dot(x.T,x)
array([[  1.,   0.,   0.,   0.],
       [  0.,   1.,   0.,   0.],
       [  0.,   0.,  35.,  44.],
       [  0.,   0.,  44.,  56.]])

>>> np.linalg.inv(np.dot(x.T,x))
array([[ 1.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  1.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  2.33333333, -1.83333333],
       [ 0.        ,  0.        , -1.83333333,  1.45833333]])