# How to extract the list of all matrices from a Block Diagonal Matrix?

I am working on a function to do the opposite operation of How to form a block-diagonal Matrix from a list of matrices?

Here is my current algorithm: For i starting at 1, scan the elements of index (i,1...i-1) and (1...i-1,i), if they are all zero, take the top left $(i-1)(i-1)$ matrix out and return the rest matrix from (i,i) as (1,1).

extractBDM[mat_?MatrixQ] :=
Flatten[Last@Reap@NestWhile[
Module[{pos},
pos = FirstPosition[
Table[Plus @@ (#[[i + 1, ;; i]] ~
Join~#[[;; i, i + 1]])^2, {i, 1, Length@# - 1}], 0];
If[MissingQ[pos], Sow[#]; {}, pos = pos[[1]];
Sow[#[[;; pos, ;; pos]]]; #[[pos + 1 ;;, pos + 1 ;;]]]] &,
mat, # != {} &], 1]


However this algorithm does not support the existence of non-square matrices on the diagonal, and the Table inside is calculating the whole matrix, while this is unnecessary. Are there better ways to solve this problem?

On uniqueness: In the result matrices, for every diagonal element, say it is index is {i,i}, there should be at least one non-zero element among {1,i}, {2,i}, ...,{i-1,i} and {i,1}, {i,2}, ...,{i,i-1}; also the diagonal elements are all non-zero.

• Just a toy example which is assuming there are no 0 in blocks: SplitBy[m, Unitize] /. (0 -> Nothing) // Map@MatrixForm
– Kuba
Jul 11, 2016 at 8:13
• +1 for an interesting question, but you realize that allowing zeros in the blocks, noted by Kuba, means some reversals are not unique. How would you propose those are to be handled?
– ciao
Jul 11, 2016 at 8:53
• @ciao I added some clarifications, is there any ambiguity now? Jul 11, 2016 at 9:01
• How would you want $$\begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix}$$ to be handled? Jul 11, 2016 at 9:04
• @J.M. I have added that the diagonal elements are non-zero. But This problem is not a fake one, it is from a real math problem. The matrix is a result from Jordan Decomposition, and even if there are sometimes zeros in a block, the 1 is likely to stick with other non-zero numbers. I know you are an expert in mathematics, I would appreciate it if you can help me describe this question clearly. Jul 11, 2016 at 9:09

extractBDM[mat_?MatrixQ] :=
Pick[mat, #, 1] & /@
Values[ComponentMeasurements[
Image[MorphologicalTransform[Unitize[mat], "BoundingBoxes"]],
"Mask", CornerNeighbors -> False]] /. {} -> Nothing


# Explanation

I just treat the values as an image pixel.So our purpose is finding the connected component.Supose your list mat is

MatrixForm[
mat = {{10, 2, 0, 0, 0}, {0, 5, 0, 0, 0}, {0, 0, 7, 3, 0}, {0, 0, 0,
4, 0}, {0, 0, 0, 0, 6}}]


$\left( \begin{array}{ccccc} 10 & 2 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 & 0 \\ 0 & 0 & 7 & 3 & 0 \\ 0 & 0 & 0 & 4 & 0 \\ 0 & 0 & 0 & 0 & 6 \\ \end{array} \right)$

Well these pixel of bounding boxes is what we want

Image[MorphologicalTransform[Unitize[mat], "BoundingBoxes"]]


Then we use Pick to select every component's pixel in mat

# Usage

MatrixForm[mat]


$\left( \begin{array}{ccccc} 10 & 2 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 & 0 \\ 0 & 0 & 7 & 3 & 0 \\ 0 & 0 & 0 & 4 & 0 \\ 0 & 0 & 0 & 0 & 6 \\ \end{array} \right)$

extractBDM[mat]


{{{10,2},{0,5}},{{7,3},{0,4}},{{6}}}

MatrixForm[
m = {{5, 3, 0, 0, 0, 0, 0, 0}, {4, 2, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 9,
9, 2, 0, 0}, {0, 0, 4, 2, 3, 6, 0, 0}, {0, 0, 3, 7, 5, 1, 0,
0}, {0, 0, 0, 0, 0, 0, 3, 1}, {0, 0, 0, 0, 0, 0, 1, 7}}]


$\left( \begin{array}{cccccccc} 5 & 3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 4 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 9 & 9 & 2 & 0 & 0 \\ 0 & 0 & 4 & 2 & 3 & 6 & 0 & 0 \\ 0 & 0 & 3 & 7 & 5 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 3 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 7 \\ \end{array} \right)$

MatrixForm /@ extractBDM[m]


• +1 I just knew there had to be a graphical solution Jul 12, 2016 at 7:56
• @LLlAMnYP Thanks.Just strike to my mind. :)
– yode
Jul 12, 2016 at 8:14
• "+1", but this method doesn't work if the matrix is huge and the blocks are a bit sparse.
– Ali
Jan 14, 2018 at 1:17

I am not sure if this is what you want (it would be helpful if you can give an example in your question). What I am doing here is to get the block matrix at any given position.

First I create a general matrix.

a = RandomInteger[{1, 9}, {2, 2}];
b = RandomInteger[{1, 9}, {3, 4}];
c = RandomInteger[{1, 9}, {2, 2}];
m = SparseArray[{Band[{1, 1}] -> a, Band[{3, 3}] -> b, Band[{6, 7}] -> c}] // Normal;
m = ReplacePart[m, {3, 3} -> 0];
MatrixForm[m]


$\begin{pmatrix} 5 & 3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 4 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 9 & 9 & 2 & 0 & 0 \\ 0 & 0 & 4 & 2 & 3 & 6 & 0 & 0 \\ 0 & 0 & 3 & 7 & 5 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 3 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 7 \\ \end{pmatrix}$

Now I want to get the block at {3,3}

{n1, n2} = {3, 3};

r1 = Flatten[SparseArray[m[[n1, All]]]["NonzeroPositions"]];
c1 = Flatten[SparseArray[m[[All, n2]]]["NonzeroPositions"]];
r2 = Flatten[SparseArray[m[[#, All]]]["NonzeroPositions"] & /@ c1];
c2 = Flatten[SparseArray[m[[All, #]]]["NonzeroPositions"] & /@ r1];
{b21, b22} = {Min[#], Max[#]} &@c2;
{b11, b12} = {Min[#], Max[#]} &@r2;
m1 = m[[b21 ;; b22, b11 ;; b12]];

MatrixForm[m1]


$\begin{pmatrix} 0 & 9 & 9 & 2 \\ 4 & 2 & 3 & 6 \\ 3 & 7 & 5 & 1 \\ \end{pmatrix}$

For the first block take {n1, n2} = {1, 1}.

• Thanks, but your function does not work properly when the last elements of the first row/column are zero, eg. the first row is 9,9,0. Jul 11, 2016 at 8:32
• @happyfish, now it will work. I modified it for any diagonal block. Jul 11, 2016 at 9:04

I'll give two solutions here:

1. a not-that elegant solution which can apply to multiple sizes but all blocks MUST be Fully filled! (Just the limitations of Morphological methods)

Function[{n}, Pick[mat, #, n] /. {} -> Nothing] /@ Range@Max@# &@
MorphologicalComponents[mat, Method -> "BoundingBox",
CornerNeighbors -> False]


This use mainly the MorphologicalComponents function and find out the blocks in a Morphological way.

1. If, as I suppose, you're dealing with Jordan Blocks :) Then this will work.

mat[[#, #]] & /@ Thread[Prepend[#, 1] ;; Append[#, Length@mat]] &@
Flatten@Position[Diagonal[mat, 1], 0]


This use the property of Jordan blocks. Check MathWorld and find it out~

# UPDATE for Solution 1

The limitation of MorphologicalComponents can be overcome with some easy manipulation. Changing it to this function and it will always do a great job~

newMorphologicalComponents[mat_] :=
NestWhile[
MorphologicalComponents[#, CornerNeighbors -> False,
Method -> "BoundingBox"] &, mat, UnsameQ, 2]


As you can see, this will do multiple times of Morphological components till all overlapping components merge into one huge cluster which is just what we need.

• the first code is, in some aspects, a simplified and more direct solution of @yode 's solution.
– Wjx
Jul 12, 2016 at 14:11
• Poor composing but actually have better efficiency than me. :)
– yode
Jul 13, 2016 at 17:43
• "than my solution" or "than mine"~ :P
– Wjx
Jul 14, 2016 at 0:42
• Aha..Thanks for pointing that.
– yode
Jul 14, 2016 at 1:23
• Updated, this time MorphologicalComponents' drawback has been overcame~@yode
– Wjx
Jul 15, 2016 at 9:46

You can treat your input matrix as the adjacency matrix of a graph, and then each block corresponds to a different connected component in the graph. You can conveniently extract these using the function ConnectedComponents[], so a one line answer to your question is then

Table[m[[comp, comp]], {comp, ConnectedComponents@AdjacencyGraph@Unitize@m}]


or, more succinctly,

m[[#, #]] & /@ ConnectedComponents@AdjacencyGraph@Unitize@m

• Note there are some thing should be tunning still.
– yode
Jul 22, 2017 at 10:10
• It is unclear to me how to define a (unique and reasonable) decomposition when blocks are non-square. For example, the matrix { {1, 1, 1, 0, 0}, {1, 1, 1, 0, 0}, {0, 0, 1, 1, 1}, {0, 0, 0, 1, 1}, {0, 0, 0, 1, 1} } may be decomposed in a number of different ways, in two or three blocks, etc Nov 30, 2017 at 5:57
• Another thing to be aware of is the case where connected components are between nods which aren't "close" to each other. For example, if you only have connections form node 5 to 8, and only between 6 and 7, then this algorithm will output a block for 5 and 8 and another one for 6,7. To fix it just take the connected components and add the missing numbers, {5,8} --> {5,6,7,8}. Or reorder the basis and avoid all this. Depends on how flexible the problem is. Sep 20, 2020 at 10:11