In 1858, in A Memoir on the Theory of Matrices, Arthur Cayley systematically treated the matrix as a complete algebraic entity—not a collection of coefficients taken from a system of equations, but a “single object” that can be added, multiplied, and inverted. This shift in viewpoint leads naturally to a question: if this object is overwhelmingly large, can we cut it into pieces and perform algebraic operations directly on those pieces, just as we would with ordinary numbers?
Over the following century, this question spread from the halls of abstract algebra to the computational work of science and engineering, and in the end it profoundly reshaped the underlying architecture of modern computing.
The Schur complement shows the power of block operations. Let the block matrix with compatible sizes be
If is an invertible square matrix, block Gaussian elimination gives . Here is called the Schur complement. When the sizes are compatible and the blocks that must be inverted are invertible, submatrices can serve as units of computation; but matrix multiplication is in general not commutative, so the order of the factors cannot be changed at will. This structure has important applications in multivariate Gaussian distributions, control theory, and convex optimization.
As problems in science and engineering grew in scale, “blocking” also became an important way to organize computation: divide a large matrix into manageable sub-blocks, schedule the movement of data, distribute the work, and then combine the results according to compatible rules of operation. Today’s distributed computing and CPU cache blocking still exploit this idea. The concrete benefits depend on the data layout, the hardware, and the algorithm; partitioning into blocks does not by itself guarantee fewer arithmetic operations.
In 1969, Volker Strassen showed that partitioning into blocks can also change the number of operations. The standard algorithm needs scalar multiplications to compute the product of two matrices, but this is not a proven lower bound. He constructed a formula for block multiplication that needs only 7 block multiplications instead of the 8 required by direct expansion.
Precisely because the objects being operated on inside the blocks are still matrices, the technique can be applied recursively down to scalars or to suitably small blocks, lowering the computational complexity of matrix multiplication directly from to . Strassen’s paper bore a strikingly provocative title: Gaussian Elimination is not Optimal. This was not merely a matter of saving a few multiplications; it changed how people understood the operation count of the standard algorithm and set off a “gold rush” around the matrix multiplication exponent that is still being fought out at the research frontier today.
The story of block matrices teaches us that structure is power. When the dimension of a problem is hopelessly high, seeing the block geometry inside it is like holding the key that cuts the problem down to size. Partitioning into blocks is not a passive compromise but a profound way of thinking that makes complex large-scale systems modular and parallel—from the stiffness assembly of finite element analysis and the partial trace of density matrices in quantum mechanics to the tensor blocking of the self-attention mechanism in today’s deep-learning Transformer architectures, this key keeps turning in the deepest layers of the computational world.
Chapter Structure and Learning Objectives¶
Cayley’s idea of “the matrix as an object” is the prerequisite for block thinking, the Schur complement is an important application of it, and Strassen’s algorithm is its most counterintuitive demonstration. The task of this chapter is to place this way of thinking on a rigorous algebraic foundation.
§5.1 first establishes the language: the definition of block matrices, together with addition, scalar multiplication, and transposition. The rules for these basic operations are entirely analogous to those for ordinary matrices, but they call for care with one extra condition, the “compatibility of block sizes.”
§5.2 is the core of the chapter: block multiplication. Two block matrices are multiplied by a rule that has exactly the same form as ordinary matrix multiplication—simply replace “entries” by “submatrix blocks”—provided that the dimensions of adjacent blocks match. Strassen’s algorithm was born within the framework of this section; its formula with 7 multiplications looks like magic, yet careful checking shows it to be exactly right.
§5.3 discusses inversion of block matrices. For a block structure, the Schur complement provides an elegant inversion formula and at the same time reveals the precise relationship between “the large matrix is invertible” and “certain sub-blocks are invertible.”
§5.4 turns to advanced topics: the tensor product (Kronecker product) and the partial trace. The tensor product “weaves” two matrices into a larger matrix and is the mathematical foundation both of composite systems in quantum mechanics and of the multi-head attention mechanism in deep learning; the partial trace is a reduction operation that eliminates the degrees of freedom of one subsystem, is in general not invertible, and is used to “extract” information about a subsystem from a composite system.
Once you have read this chapter, your fear of “large matrices” should be somewhat smaller—because you will know that a large matrix is nothing more than a few small matrices put together according to a pattern.
5.1 Definition and Basic Operations of Block Matrices¶
A block matrix is a large matrix divided into several smaller matrix blocks, which can be regarded as the “sub-entries” of the matrix. Partitioning in this way turns operations on a large matrix into operations on the individual small matrices, thereby simplifying the computation.
Below we introduce three basic operations on block matrices.
5.1.1 Addition of Block Matrices¶
If two block matrices have the same block structure, their sum is defined by adding corresponding blocks. Let
Then
5.1.2 Scalar Multiplication of Block Matrices¶
To perform scalar multiplication on a block matrix, simply multiply each block by the scalar. That is, if is a scalar,
5.1.3 Transpose of Block Matrices¶
Taking the transpose of a block matrix requires both transposing each block and swapping the positions of the blocks. Writing for the transpose of , we have
5.1.4 ◆4D Tensor Representation and Indexing of Block Matrices¶
Besides being viewed as a two-dimensional structure made up of small matrix blocks, a block matrix can also be regarded as a 4D tensor; this viewpoint is important in modern machine learning and quantum mechanics. With the 4D tensor representation we can handle batched block matrix operations more naturally and take advantage of the tensor operations offered by modern computing frameworks.
For a matrix partitioned into a block structure, we have the following correspondence:
2D matrix index: , where
4D tensor index: , where
(block-row index)
(block-column index)
(within-block row index)
(within-block column index)
Summary.
Representing a block matrix as a 4D tensor gives us a powerful tool, especially in modern computing environments that call for large numbers of block-level operations. Understanding this representation and its indexing rules helps us design and implement block matrix algorithms more effectively.
5.1.5 Python: Block Matrix Conversions with NumPy¶
import numpy as np
import matplotlib.pyplot as plt
print("=" * 60)
print("5.1 Building block matrices - using np.block, np.stack, np.hstack")
print("=" * 60)
# ========== 1. Defining the basic blocks ==========
print("\n1. Define the basic matrix blocks")
print("-" * 30)
# Define four 2×2 submatrix blocks
A00 = np.array([[1, 2],
[3, 4]])
A01 = np.array([[5, 6],
[7, 8]])
A10 = np.array([[9, 10],
[11, 12]])
A11 = np.array([[13, 14],
[15, 16]])
print(f"A00 (upper-left block):\n{A00}")
print(f"\nA01 (upper-right block):\n{A01}")
print(f"\nA10 (lower-left block):\n{A10}")
print(f"\nA11 (lower-right block):\n{A11}")
# ========== 2. Assembling block matrices with np.block ==========
print("\n\n2. Assemble a block matrix with np.block")
print("-" * 40)
# np.block is the most intuitive way to build a block matrix
block_matrix = np.block([[A00, A01],
[A10, A11]])
print("Using np.block([[A00, A01], [A10, A11]]):")
print(block_matrix)
print(f"Shape: {block_matrix.shape}")
# A more complex block structure
print("\nComplex block structure - zero and identity matrices:")
I2 = np.eye(2) # 2×2 identity matrix
O2 = np.zeros((2, 2)) # 2×2 zero matrix
complex_block = np.block([[A00, O2, A01],
[I2, A10, O2],
[A11, I2, A00]])
print(f"Shape: {complex_block.shape}")
print(complex_block)
# Combining blocks of different sizes
print("\nCombining blocks of different sizes:")
B1 = np.ones((3, 2)) # 3×2 all-ones matrix
B2 = np.eye(3) # 3×3 identity matrix
B3 = np.zeros((2, 2)) # 2×2 zero matrix
B4 = np.full((2, 3), 5) # 2×3 matrix filled with 5
mixed_block = np.block([[B2, B1],
[B4, B3]])
print(f"Shape of the mixed-size block matrix: {mixed_block.shape}")
print(mixed_block)
# ========== 3. Horizontal stacking with np.hstack ==========
print("\n\n3. Horizontal stacking with np.hstack")
print("-" * 35)
# First assemble the top and bottom rows separately
top_row = np.hstack([A00, A01])
bottom_row = np.hstack([A10, A11])
print("Top row np.hstack([A00, A01]):")
print(top_row)
print(f"Shape: {top_row.shape}")
print("\nBottom row np.hstack([A10, A11]):")
print(bottom_row)
print(f"Shape: {bottom_row.shape}")
# Then stack vertically
hstack_result = np.vstack([top_row, bottom_row])
print("\nFinal result np.vstack([top_row, bottom_row]):")
print(hstack_result)
print(f"Shape: {hstack_result.shape}")
# Check agreement with the np.block result
print(f"\nAgrees with the np.block result: {np.array_equal(block_matrix, hstack_result)}")
# Building in several steps
print("\nBuilding a large block matrix in several steps:")
# Create more blocks
C11 = np.full((2, 2), 17)
C12 = np.full((2, 2), 18)
C21 = np.full((2, 2), 19)
C22 = np.full((2, 2), 20)
# Step 1: combine horizontally
row1 = np.hstack([A00, A01, C11, C12])
row2 = np.hstack([A10, A11, C21, C22])
print(f"Shape of block row 0: {row1.shape}")
print(f"Shape of block row 1: {row2.shape}")
# Step 2: combine vertically
large_matrix = np.vstack([row1, row2])
print(f"Shape of the large matrix: {large_matrix.shape}")
# ========== 4. Using np.stack to help build a 2D block matrix ==========
print("\n\n4. Using np.stack to help build a 2D block matrix")
print("-" * 45)
# Method 1: organize the rows with np.stack, then reshape to 2D
print("Method 1: stack into rows first, then reshape")
top_blocks = [A00, A01]
bottom_blocks = [A10, A11]
# Stack along the row direction
top_row_stack = np.stack(top_blocks, axis=1) # 2×2×2 -> reshaped to 2×4
bottom_row_stack = np.stack(bottom_blocks, axis=1) # 2×2×2 -> reshaped to 2×4
print(f"Shape of the stacked top row: {top_row_stack.shape}")
print(f"Shape of the stacked bottom row: {bottom_row_stack.shape}")
# Reshape to 2D
top_row_2d = top_row_stack.reshape(2, 4)
bottom_row_2d = bottom_row_stack.reshape(2, 4)
print("Top row reshaped to 2D:")
print(top_row_2d)
print("Bottom row reshaped to 2D:")
print(bottom_row_2d)
# Final vertical combination
stack_result = np.vstack([top_row_2d, bottom_row_2d])
print("\nFinal 2D block matrix:")
print(stack_result)
print(f"Shape: {stack_result.shape}")
# Method 2: organize everything with np.stack, then reshape
print("\nMethod 2: stack everything, then reshape")
# Arrange in the logical order of the block matrix
blocks_2x2 = np.array([[A00, A01],
[A10, A11]]) # shape: (2, 2, 2, 2)
print(f"Shape of the 4D block structure: {blocks_2x2.shape}")
# Reshape into a 2D block matrix
stack_method2 = np.concatenate([
np.concatenate([blocks_2x2[0, 0], blocks_2x2[0, 1]], axis=1),
np.concatenate([blocks_2x2[1, 0], blocks_2x2[1, 1]], axis=1)
], axis=0)
print("Reshaped 2D matrix:")
print(stack_method2)
# Method 3: handle different kinds of blocks with np.stack
print("\nMethod 3: handling blocks of mixed kinds")
# Create different kinds of 2×2 blocks
zero_block = np.zeros((2, 2))
ones_block = np.ones((2, 2))
eye_block = np.eye(2)
random_block = np.random.randint(1, 5, (2, 2))
print("Zero block:")
print(zero_block)
print("All-ones block:")
print(ones_block)
print("Identity block:")
print(eye_block)
print("Random block:")
print(random_block)
# Organize these different blocks with stack
special_blocks = [zero_block, ones_block, eye_block, random_block]
# Create a 2×2 block structure
mixed_top = np.hstack([special_blocks[0], special_blocks[1]])
mixed_bottom = np.hstack([special_blocks[2], special_blocks[3]])
mixed_result = np.vstack([mixed_top, mixed_bottom])
print("\nBlock matrix with mixed kinds of blocks:")
print(mixed_result)
print(f"Shape: {mixed_result.shape}")
# Check that all methods give consistent results
print(f"\nMethod 1 equals the original block_matrix: {np.array_equal(stack_result, block_matrix)}")
print(f"Method 2 equals the original block_matrix: {np.array_equal(stack_method2, block_matrix)}")
# ========== 5. Converting a 2D tensor into a 4D tensor ==========
print("\n\n5. Converting a 2D tensor into a 4D tensor")
print("-" * 35)
print("Original 2D block matrix:")
print(f"Shape: {block_matrix.shape}")
print(block_matrix)
# Method 1: convert the block matrix into a 4D tensor with reshape
# Assume the original matrix is 4×4 and consists of 2×2 blocks; convert it to (2,2,2,2)
tensor_4d_method1 = block_matrix.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3)
for i in range(2):
for j in range(2):
assert np.array_equal(tensor_4d_method1[i, j], block_matrix[2*i:2*i+2, 2*j:2*j+2])
print(f"\nMethod 1 - reshape: {block_matrix.shape} -> {tensor_4d_method1.shape}")
print("Block [0,0,:,:] of the 4D tensor:")
print(tensor_4d_method1[0, 0])
print("Block [0,1,:,:] of the 4D tensor:")
print(tensor_4d_method1[0, 1])
# Method 2: regroup manually into a 4D tensor, keeping the meaning of the block structure
print(f"\nMethod 2 - regrouping the block structure manually:")
# Extract each 2×2 block
block_00 = block_matrix[0:2, 0:2] # A00
block_01 = block_matrix[0:2, 2:4] # A01
block_10 = block_matrix[2:4, 0:2] # A10
block_11 = block_matrix[2:4, 2:4] # A11
# Combine into a 4D tensor (block row, block column, within-block row, within-block column)
tensor_4d_method2 = np.array([[block_00, block_01],
[block_10, block_11]])
print(f"Shape of the 4D tensor: {tensor_4d_method2.shape}")
print("Block [0,1] (corresponds to A01):")
print(tensor_4d_method2[0, 1])
print("Block [1,1] (corresponds to A11):")
print(tensor_4d_method2[1, 1])
# Method 3: rebuild the block structure by slicing
print(f"\nMethod 3 - building the 4D tensor by slicing:")
# A more flexible slicing approach, for equal-sized blocks that divide the matrix evenly
def create_4d_tensor_from_blocks(matrix, block_rows, block_cols):
"""
Convert a 2D block matrix into a 4D tensor
matrix: the input 2D matrix
block_rows, block_cols: the number of rows and the number of columns of each block
"""
matrix = np.asarray(matrix)
if matrix.ndim != 2:
raise ValueError("matrix must be a two-dimensional matrix")
if any(isinstance(x, (bool, np.bool_)) or not isinstance(x, (int, np.integer)) or x <= 0 for x in (block_rows, block_cols)):
raise ValueError("the block sizes must be positive integers")
rows, cols = matrix.shape
if rows % block_rows or cols % block_cols:
raise ValueError("both dimensions of the matrix must be divisible by the block sizes")
n_block_rows = rows // block_rows
n_block_cols = cols // block_cols
# Create the 4D tensor
tensor_4d = np.zeros((n_block_rows, n_block_cols, block_rows, block_cols), dtype=matrix.dtype)
for i in range(n_block_rows):
for j in range(n_block_cols):
start_row = i * block_rows
end_row = (i + 1) * block_rows
start_col = j * block_cols
end_col = (j + 1) * block_cols
tensor_4d[i, j] = matrix[start_row:end_row, start_col:end_col]
return tensor_4d
tensor_4d_method3 = create_4d_tensor_from_blocks(block_matrix, 2, 2)
print(f"By slicing: {block_matrix.shape} -> {tensor_4d_method3.shape}")
print("Block [0,0] (corresponds to A00):")
print(tensor_4d_method3[0, 0])
print("Block [1,1] (corresponds to A11):")
print(tensor_4d_method3[1, 1])
# Method 4: partition into blocks with np.split
print(f"\nMethod 4 - partitioning with np.split:")
# First split by rows
row_splits = np.split(block_matrix, 2, axis=0) # split into 2 rows
print(f"Splitting by rows gives {len(row_splits)} submatrices, each of shape: {row_splits[0].shape}")
# Then split each row by columns
blocks_list = []
for row_block in row_splits:
col_splits = np.split(row_block, 2, axis=1) # split each row into 2 columns
blocks_list.append(col_splits)
# Convert to a 4D tensor
tensor_4d_method4 = np.array(blocks_list)
print(f"Using split: {block_matrix.shape} -> {tensor_4d_method4.shape}")
print("Block [0,0] (corresponds to A00):")
print(tensor_4d_method4[0, 0])
print("Block [1,1] (corresponds to A11):")
print(tensor_4d_method4[1, 1])
# ========== 6. Converting a 4D tensor back into a 2D tensor ==========
print("\n\n6. Converting a 4D tensor back into a 2D tensor")
print("-" * 30)
print("Now we convert the 4D tensor back into a 2D block matrix")
# Use the 4D tensor generated above (taking method2 as an example)
print(f"Shape of the original 4D tensor: {tensor_4d_method2.shape}")
print("Contents of the 4D tensor:")
print("Block [0,0]:", tensor_4d_method2[0,0])
print("Block [0,1]:", tensor_4d_method2[0,1])
print("Block [1,0]:", tensor_4d_method2[1,0])
print("Block [1,1]:", tensor_4d_method2[1,1])
# Method 1: restore the axis order with transpose, then reshape
print(f"\nMethod 1 - using reshape:")
# Reshape from (2,2,2,2) to (4,4)
reconstructed_1 = tensor_4d_method1.transpose(0, 2, 1, 3).reshape(4, 4)
print(f"Shape conversion: {tensor_4d_method1.shape} -> {reconstructed_1.shape}")
print("Reconstructed 2D matrix:")
print(reconstructed_1)
# Method 2: using swapaxes + reshape
print(f"\nMethod 2 - using swapaxes + reshape:")
# Swap axes first, then reshape
swapped = tensor_4d_method2.swapaxes(1, 2) # (2,2,2,2) -> (2,2,2,2), but rearranged
print(f"Shape after swapping axes: {swapped.shape}")
reconstructed_2 = swapped.reshape(4, 4)
print("Reconstructed 2D matrix:")
print(reconstructed_2)
# Method 3: regroup manually (the clearest method)
print(f"\nMethod 3 - regrouping the blocks manually:")
# Extract each block
block_00 = tensor_4d_method3[0, 0] # upper left
block_01 = tensor_4d_method3[0, 1] # upper right
block_10 = tensor_4d_method3[1, 0] # lower left
block_11 = tensor_4d_method3[1, 1] # lower right
# Regroup with hstack and vstack
top_row = np.hstack([block_00, block_01])
bottom_row = np.hstack([block_10, block_11])
reconstructed_3 = np.vstack([top_row, bottom_row])
print("Reconstructed 2D matrix:")
print(reconstructed_3)
# Method 4: using np.block (the most concise)
print(f"\nMethod 4 - using np.block:")
reconstructed_4 = np.block([[tensor_4d_method4[0,0], tensor_4d_method4[0,1]],
[tensor_4d_method4[1,0], tensor_4d_method4[1,1]]])
print("Reconstructed 2D matrix:")
print(reconstructed_4)
# Method 5: using concatenate
print(f"\nMethod 5 - using concatenate:")
# First join the blocks of each row along axis 3 (within-block columns)
row_0 = np.concatenate([tensor_4d_method2[0,0], tensor_4d_method2[0,1]], axis=1)
row_1 = np.concatenate([tensor_4d_method2[1,0], tensor_4d_method2[1,1]], axis=1)
# Then join all the rows along axis 2 (within-block rows)
reconstructed_5 = np.concatenate([row_0, row_1], axis=0)
print("Reconstructed 2D matrix:")
print(reconstructed_5)
# General function: 4D tensor to 2D matrix
def tensor_4d_to_2d(tensor_4d):
"""
Convert a 4D tensor (block row, block column, within-block row, within-block column) into a 2D matrix
"""
n_block_rows, n_block_cols, block_rows, block_cols = tensor_4d.shape
# Approach: use reshape together with an axis permutation
# First permute to (block row, within-block row, block column, within-block column)
reshaped = tensor_4d.transpose(0, 2, 1, 3)
# Then reshape to 2D
result = reshaped.reshape(n_block_rows * block_rows, n_block_cols * block_cols)
return result
print(f"\nTesting the general function:")
reconstructed_general = tensor_4d_to_2d(tensor_4d_method2)
print("2D matrix reconstructed by the general function:")
print(reconstructed_general)
# Check the correctness of all methods
print(f"\nVerification results:")
print(f"Original matrix:")
print(block_matrix)
all_methods = [reconstructed_1, reconstructed_3, reconstructed_4, reconstructed_5, reconstructed_general]
method_names = ["reshape", "manual regrouping", "np.block", "concatenate", "general function"]
for i, (method_result, name) in enumerate(zip(all_methods, method_names)):
is_equal = np.array_equal(block_matrix, method_result)
print(f"{name}: equals the original matrix = {is_equal}")
# Special case: handling 4D tensors of other sizes
print(f"\nHandling 4D tensors of other sizes:")
# Create an example 4D tensor with 3×3 blocks
large_blocks = []
for i in range(3):
row_blocks = []
for j in range(3):
# each block is 2×2
block = np.full((2, 2), i*3 + j + 1)
row_blocks.append(block)
large_blocks.append(row_blocks)
large_4d = np.array(large_blocks)
print(f"Shape of the large 4D tensor: {large_4d.shape}")
# Convert to a 6×6 2D matrix
large_2d = tensor_4d_to_2d(large_4d)
print(f"Shape of the converted 2D matrix: {large_2d.shape}")
print("Converted 2D matrix:")
print(large_2d)
print("\n" + "=" * 60)
print("Summary of methods for building and converting block matrices:")
print("1. np.block: the most intuitive way to assemble blocks")
print("2. np.hstack/vstack: horizontal/vertical stacking, suited to step-by-step construction")
print("3. np.stack: stacking along a new axis, creating a structure for batch processing")
print("4. 2D→4D conversion: reshape, manual regrouping, slicing")
print("=" * 60)5.2 Multiplication of Block Matrices¶
Matrix multiplication admits four complementary geometric views, each revealing a different computational structure. Understanding these four views is the key both to understanding “why” general matrix multiplication is defined as it is and to extending multiplication to the language of blocks.
5.2.1 View One: Inner Products—Entry-by-Entry Computation¶
where denotes row of and denotes column of . Each entry of the result matrix is the inner product of a row of the left matrix with a column of the right matrix.
5.2.2 View Two: Linear Combinations of the Column Vectors of ¶
Each column of is a linear combination of the column vectors of , with coefficients supplied by the corresponding column of :
In other words, column of is the linear combination of the columns of whose coefficients are the components of .
Batched generalization (block view). If the columns of are partitioned into several blocks , then
Each block can be computed independently in parallel—this is precisely the algebraic basis of batched matrix operations on GPUs.
5.2.3 View Three: Linear Combinations of the Row Vectors of ¶
Symmetrically, each row of is a linear combination of the row vectors of , with coefficients supplied by the corresponding row of :
In other words, row of is the linear combination of the rows of whose coefficients are the components of .
5.2.4 Block Matrix Multiplication¶
In the language of block matrices, the three views above are unified by a single formula.
5.2.5 View Four: Outer-Product Expansion—Sums of Rank-One Matrices¶
is called the outer product of and ; when both vectors are nonzero, it is a rank-one matrix.
The block multiplication theorem has a particularly intuitive geometric interpretation: when the blocks degenerate into single columns of the left factor and single rows of the right factor, each term is a matrix of rank at most 1; when both vectors are nonzero its rank is 1, and the whole product is the sum of these rank-one matrices:
Generalizing to the language of blocks, each term is a “block-product contribution” whose rank need not be 1; summing them recovers the full product—this is precisely the geometric essence of the theorem.
5.2.6 ◆ Strassen’s Algorithm: The Power of Block Multiplication¶
A striking application of the block multiplication theorem is Strassen’s algorithm (Strassen, 1969). Intuitively, multiplying two matrices takes 8 scalar multiplications; Strassen discovered that 7 suffice, by constructing the following 7 intermediate quantities:
and then recombining them using additions:
Applying this technique recursively to matrices lowers the computational complexity from to .
5.2.7 Python: Block Matrix Multiplication with NumPy¶
import numpy as np
print("Comparing block matrix multiplication with 4D tensors")
print("=" * 40)
# Create block matrices in 4D tensor form
# A: (2, 3, 4, 4) - 2x3 blocks, each a 4x4 matrix
# B: (3, 2, 4, 4) - 3x2 blocks, each a 4x4 matrix
np.random.seed(42)
A_4d = np.random.randn(2, 3, 4, 4)
B_4d = np.random.randn(3, 2, 4, 4)
print(f"Shape of A as a 4D tensor: {A_4d.shape}")
print(f"Shape of B as a 4D tensor: {B_4d.shape}")
print("\n" + "-" * 40)
print("Method 1: convert back to 2D matrices, then multiply as usual")
print("-" * 40)
def tensor_4d_to_2d(tensor_4d):
"""Convert a 4D block tensor into a 2D matrix"""
num_block_rows, num_block_cols, block_height, block_width = tensor_4d.shape
# Compute the size of the final matrix
total_height = num_block_rows * block_height
total_width = num_block_cols * block_width
# Initialize the result matrix
matrix_2d = np.zeros((total_height, total_width), dtype=tensor_4d.dtype)
# Fill in each block
for i in range(num_block_rows):
for j in range(num_block_cols):
row_start = i * block_height
row_end = (i + 1) * block_height
col_start = j * block_width
col_end = (j + 1) * block_width
matrix_2d[row_start:row_end, col_start:col_end] = tensor_4d[i, j]
return matrix_2d
# Convert the 4D tensors into 2D matrices
A_2d = tensor_4d_to_2d(A_4d)
B_2d = tensor_4d_to_2d(B_4d)
print(f"Shape of A as a 2D matrix: {A_2d.shape}")
print(f"Shape of B as a 2D matrix: {B_2d.shape}")
# Ordinary matrix multiplication
C_2d_method1 = A_2d @ B_2d
print(f"Shape of the method 1 result: {C_2d_method1.shape}")
print("\n" + "-" * 40)
print("Method 2: multiply the 4D tensors with einsum, then convert back")
print("-" * 40)
# Block matrix multiplication with einsum
# 'ikmn,klno->ilmo': k and n are the summation indices
# i,l: block-row and block-column indices of the result; m,o: within-block row and column indices
C_4d = np.einsum('ikmn,klno->ilmo', A_4d, B_4d)
print(f"Shape of the 4D einsum result: {C_4d.shape}")
# Convert the result into a 2D matrix
C_2d_method2 = tensor_4d_to_2d(C_4d)
print(f"Shape of the method 2 result: {C_2d_method2.shape}")
print("\n" + "-" * 40)
print("Comparing the results")
print("-" * 40)
# Compare the results of the two methods
difference = np.abs(C_2d_method1 - C_2d_method2)
max_diff = np.max(difference)
mean_diff = np.mean(difference)
print(f"The two methods agree: {np.allclose(C_2d_method1, C_2d_method2)}")
print(f"Maximum difference: {max_diff:.2e}")
print(f"Mean difference: {mean_diff:.2e}")
print("\n" + "-" * 40)
print("Detailed check: computing one block by hand")
print("-" * 40)
# Compute block (0,0) of the result by hand as a check
# C[0,0] = A[0,0]@B[0,0] + A[0,1]@B[1,0] + A[0,2]@B[2,0]
manual_block_00 = (A_4d[0,0] @ B_4d[0,0] +
A_4d[0,1] @ B_4d[1,0] +
A_4d[0,2] @ B_4d[2,0])
# Extract the corresponding block from the method 1 result
method1_block_00 = C_2d_method1[0:4, 0:4]
# Extract the corresponding block from the method 2 result
method2_block_00 = C_2d_method2[0:4, 0:4]
print("By hand vs method 1 vs method 2:")
print(f"By hand equals method 1: {np.allclose(manual_block_00, method1_block_00)}")
print(f"By hand equals method 2: {np.allclose(manual_block_00, method2_block_00)}")
print(f"Method 1 equals method 2: {np.allclose(method1_block_00, method2_block_00)}")5.3 Inverses of Matrices in Block Form¶
For matrices with certain special structures (such as diagonal, upper triangular, or lower triangular matrices), we can use the block structure to invert in stages and reduce the computational complexity. In the partitions considered here, every block on the diagonal is required to be a square matrix.
5.3.1 Inverses of Block Diagonal Matrices¶
These properties make computations with block diagonal matrices very efficient, because we can operate on each diagonal block separately and then combine the results.
5.3.2 The Block Formula for the Inverse and the Schur Complement¶
For a general block matrix
where is a square matrix and is a square matrix, we would like to follow the strategy of §5.3.1 and reduce the computation of to inverting smaller matrices. The difficulty is that the off-diagonal blocks and couple the blocks to one another—the way out is to first use block multiplication to “eliminate” one of the off-diagonal blocks.
Derivation: block Gaussian elimination. Assume that is invertible. Imitating the step in ordinary Gaussian elimination of “using row 0 to eliminate the entries below it,” we left-multiply by a block lower triangular matrix to eliminate the block in the lower-left corner of :
where the matrix that appears naturally in the lower-right corner,
is called the Schur complement of in . What remains after the elimination is a block upper triangular matrix, and by Proposition 1, to be established in §5.3.3 (a block triangular matrix is invertible each diagonal block is invertible), as long as is also invertible, the whole matrix is invertible.
Note the cancellation in the computation of the lower-left block—the Schur complement is precisely the “recipe” that makes all the cross terms cancel exactly. This is no coincidence; it is the trace that block Gaussian elimination leaves behind at the level of the formula.
Summary.
The Schur complement is not a formula to be memorized; it is the natural product of block Gaussian elimination. The residue left in the lower-right corner after has been eliminated measures exactly “how much invertibility retains once the coupling transmitted through has been subtracted.” Through the block formula, inverting a large matrix is reduced to inverting two smaller matrices ( and ) plus matrix multiplications—this is the embryonic form of the block elimination and LU decomposition of Chapter 6.
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)
print("=" * 60)
print("5.3.2 Schur complement block inversion formula - numerical check")
print("=" * 60)
# The four sub-blocks from the example
A = np.array([[2., 1.], [1., 1.]])
B = np.eye(2)
C = np.array([[1., 2.], [0., 1.]])
D = np.array([[3., 1.], [1., 2.]])
M = np.block([[A, B], [C, D]])
print("Block matrix M =")
print(M)
# Step 1: Schur complement S = D - C A⁻¹ B
A_inv = np.linalg.inv(A)
S = D - C @ A_inv @ B
print("\nSchur complement S = D - C A⁻¹ B =")
print(S)
# Step 2: block inversion formula
S_inv = np.linalg.inv(S)
M_inv_block = np.block([
[A_inv + A_inv @ B @ S_inv @ C @ A_inv, -A_inv @ B @ S_inv],
[-S_inv @ C @ A_inv, S_inv ]
])
print("\nM⁻¹ from the block formula =")
print(M_inv_block)
# Step 3: compare with direct inversion in NumPy
M_inv_direct = np.linalg.inv(M)
print("\nnp.linalg.inv(M), direct inversion =")
print(M_inv_direct)
print("\nThe two agree:", np.allclose(M_inv_block, M_inv_direct))
print("M @ M⁻¹ = I:", np.allclose(M @ M_inv_block, np.eye(4)))
# Step 4: Schur determinant formula det M = det A · det S
print("\ndet(M) =", np.linalg.det(M).round(4))
print("det(A) * det(S) =", (np.linalg.det(A) * np.linalg.det(S)).round(4))
5.3.3 Inverses of Upper and Lower Triangular Matrices¶
Triangular matrices are an important class of matrices with special structure, and both the computation and the properties of their inverses have distinctive features.
Basic Definitions
Basic Properties of Triangular Matrices
Triangular matrices are well behaved under multiplication: the product of two upper triangular matrices is again upper triangular, and the product of two lower triangular matrices is again lower triangular.
Inverses of Block Triangular Matrices
Continuing the block matrix inversion of the previous section, a square block matrix with triangular structure has a comparatively simple inverse, which likewise preserves the triangular structure:
Existence of Inverses of Triangular Matrices
Here is a concrete example:
Group Structure
Summary
The inversion of triangular matrices follows clear rules: invertibility is completely determined by the diagonal entries, and the inverse preserves the same triangular structure. The block form of the inversion formula provides an effective tool for handling large, complicated matrices, and the group structure lays an important foundation for matrix theory and numerical methods.
5.4 ◆Advanced Topics in Block Matrices¶
5.4.1 The Tensor Product (Kronecker Product) and Block Matrices¶
The tensor product is an important application of block matrices, especially in quantum mechanics and multilinear algebra. Chapter 4 (Definition 23) already defined the tensor product space at the level of abstract vector spaces and announced that the details of its operations would be left to this chapter. In the setting of matrix algebra, the tensor product and the Kronecker product are two names for the same concept—indeed, the space of matrices is itself isomorphic to a tensor product space, (each rank-one matrix corresponds to the pure tensor ), so the Kronecker product of this section is precisely the realization of the abstract tensor product in concrete matrix coordinates. Tensor product matrices have a natural block structure.
These properties are very useful in practical computation—especially the mixed-product property, which allows us to break complicated tensor product operations down into simpler matrix multiplications.
As the definition shows, the tensor product naturally forms a block matrix structure in which each block is the original matrix multiplied by the corresponding entry of .
Connection with the 4D tensor representation.
Note that this block structure of the tensor product corresponds exactly to the 4D tensor representation discussed in Section 5.1.4. If the tensor product is viewed as a 4D tensor, then
This corresponds to the notation of Section 5.1.4. Each block can be understood as a “slice” of the 4D tensor.
Tensor Products of Vectors¶
When we work with vectors, the tensor product can likewise be represented by a block structure. If and , then :
This shows that the tensor product of vectors can be viewed as a special kind of block matrix. In the framework of 4D tensors, it can be regarded as a degenerate case in which some axes have size 1.
Advantages of the Block Matrix Representation¶
From a computational point of view, the block representation of the tensor product has several important advantages:
Clear structure. The block representation shows directly how the tensor product is constructed and corresponds naturally to the 4D tensor indexing of Section 5.1.4
Computational convenience. The block structure can be used to design more efficient parallel algorithms, especially in combination with batched operations on 4D tensors
Memory management. For large sparse matrices, the block representation can save storage space
Numerical stability. Block operations help control the accumulation of numerical errors
Applications in Quantum Mechanics¶
In quantum mechanics, the tensor product is the basic tool for describing composite quantum systems. When two quantum systems are combined, the total state space is the tensor product of the state spaces of the subsystems. Combined with the 4D tensor representation of Section 5.1.4, we can handle these high-dimensional structures more efficiently, especially when computing the various operations involved in quantum entanglement and quantum information processing.
5.4.2 Trace and Partial Trace of Block Matrices¶
The trace and the partial trace are important when dealing with matrices that have a tensor product structure, especially in quantum information and the analysis of many-body systems. Building on the 4D tensor representation of Section 5.1.4, this section systematically introduces how these operations are computed.
5.4.2.1 The Trace of a Block Matrix¶
For tensor product matrices, the trace has a particularly simple property:
5.4.2.2 Definition and Computation of the Partial Trace¶
The partial trace is the natural generalization of the trace to tensor product structures, used to take a trace “partially.” Using the 4D tensor representation of Section 5.1.4, we can give a precise definition of the partial trace.
Computing the Partial Trace from the Block Matrix Viewpoint¶
From a practical computational standpoint, block matrices provide a more intuitive method:
Correspondence between the 4D Tensor and Block Computations¶
A Complete Worked Example¶
5.5 Chapter Summary¶
Review of the Theoretical Thread¶
This chapter started from a simple question: when a matrix is too large to handle directly, can its internal structure become the way forward for computation? The four sections unfolded in turn along a single logical thread, each answering the question left open by the one before.
§5.1 laid the foundation of the language. Although the rules for addition, scalar multiplication, and transposition of block matrices are entirely analogous to those for ordinary matrices, the extra condition of “compatibility of block sizes” runs through all of them and cannot be ignored. The 4D tensor viewpoint further gives the block indices and the within-block entry indices a single unified address, which both corresponds directly to batched computation in NumPy and lays down the conceptual interface for the Kronecker product later on.
§5.2 revealed the true power of the language of blocks. Matrix multiplication has four complementary geometric intuitions—inner products, column combinations, row combinations, and outer-product expansion—and all four converge in the block framework to one and the same formula: . This formula, which formally “changes nothing,” nevertheless enabled Strassen, with a block technique, to lower the computational complexity of matrix multiplication from to , overturning an intuition a century and a half old.
§5.3 built on multiplication to discuss strategies for inverting in blocks. From the block-by-block inversion of block diagonal matrices, to the central role played by the Schur complement in a general block structure, to the rigorous characterization of triangular matrices, for which nonzero diagonal entries are the sole condition for invertibility, all three make the same point: the structure of inversion depends on the structure of multiplication. The invertible upper triangular matrices of a fixed size form a group under multiplication, and this algebraic property directly drives the correctness argument for the LU decomposition in numerical linear algebra.
§5.4 extended block thinking to higher dimensions. The Kronecker product scales the whole of by each entry of , used as a scalar coefficient, constructing the complete state space of a composite system; when the evolution factors as , the mixed-product property allows the local actions to be computed separately, whereas general interactions admit no such factorization. The partial trace, in turn, performs a reduction, projecting the statistical information of a single subsystem out of the high-dimensional description of the composite system; in the language of 4D tensors this operation is nothing more than “summing the diagonal entries along specified axes,” echoing the viewpoint of §5.1 and closing the conceptual loop of the whole chapter.
Connections to Other Chapters¶
The logical starting point of this chapter is Chapter 4’s abstract discussion of the structure of linear spaces. Abstract linear spaces established the equivalence among vectors, linear mappings, and matrices; block matrices are a “coarse-graining” of this equivalence—binding several dimensions together into one algebraic unit, so that a matrix of matrices becomes a legitimate mathematical object. The dimension-matching condition of block multiplication is precisely the concrete expression, in the language of blocks, of the compatibility of spaces that allows linear mappings to be composed.
Chapter 6 will take up the systematic study of systems of linear equations, and §5.3 has already laid the groundwork for its most important computational tools. The Schur complement appears in statistics as the computational core of partial correlation coefficients and in control theory as the solution structure of the Riccati equation, and above all it is the algebraic essence of the block form of Gaussian elimination. The invertibility theorem and the group structure of triangular matrices correspond directly to the correctness guarantees of the forward substitution and back substitution algorithms. This connection will reappear in a clearer form when Chapter 6 discusses the LU decomposition.
The Role of This Chapter in the Book¶
The core problem that block matrices solve is how to turn a problem that is “too large” into one that “can be handled.” This problem is more universal than linear algebra itself: during World War II, human computers processed the strength matrices of aircraft wings by hand in batches; today, GPU parallel computing carries out the gradient updates of deep learning in batches; and quantum computers store the density matrices of composite systems in a tensor product structure—the same way of thinking reappears in different guises across different eras.
The shock of Strassen’s algorithm lies not only in the few multiplications it saved but even more in what it tells us: even for a “trivial” operation, once its internal structure is seen, a cleverer path may exist. This insight directly gave rise to the research area of “lower bounds on the complexity of matrix multiplication,” which remains open to this day. If you continue reading with the block viewpoint built in this chapter, you will find that the determinant expansions, eigendecompositions, and singular value decompositions of later chapters are all, in some sense, extensions of one and the same strategy: “see the structure, exploit the structure.”