Have you ever seen the "brain teaser" for children that shows a 4 x 4 grid and asks "how many squares of any size are in this grid?" To solve this problem, the reader must recognize that there are sixteen 1 x 1 squares, nine 2 x 2 squares, four 3 x 3 squares, and one 4 x 4 square. Hence there are a total of (16+9+4+1) = 30 squares of any size.

Recently a SAS statistical programmer asked a similar question about submatrices of a matrix. Specifically, the programmer wanted to form square submatrices of consecutive rows and columns.
This article shows a few tricks in the SAS/IML language that enable you to generate and compute with all
k x k submatrices of a matrix, where the submatrices are formed by using consecutive rows and columns of the original matrix. Although it is possible to consider more general submatrices, in this article, a "submatrix" always refers to one that is formed by using **consecutive rows and columns**.

### Use subscripts to form submatrices

I've previously written about several ways to extract submatrices from a matrix.
Recall the "colon operator" (formally called the *index creation operator*) generates a vector of consecutive integers. You can use the vectors in the row and column position of the subscripts to extract a submatrix. For example, the following SAS/IML statements define a 4 x 4 matrix and extract the four 3 x 3 submatrices:

A = {4 3 1 6, 2 4 3 1, 0 2 4 3, 5 0 2 4}; /* extract the four 3 x 3 submatrices */ A11 = A[1:3, 1:3]; A12 = A[1:3, 2:4]; A21 = A[2:4, 1:3]; A22 = A[2:4, 2:4]; |

Manually specifying the row and column numbers is tedious.
Let's write a SAS/IML function that will return a k x k submatrix of A by starting at the (i,j)th position of A. Call k the *order* of the submatrix.

/* Given A, return the square submatrix of order k starting at index (i,j) */ start submat(A, i, j, k); return A[ i:(i+k-1), j:(j+k-1) ]; finish; /* example : 3 x 3 submatrices of A */ A11 = submat(A, 1, 1, 3); A12 = submat(A, 1, 2, 3); A21 = submat(A, 2, 1, 3); A22 = submat(A, 2, 2, 3); /* https://blogs.sas.com/content/iml/2015/12/02/matrices-graphs-gridded-layout.html */ ods layout gridded columns=4 advance=table; print A11, A12, A21, A22; ods layout end; |

### The number of submatrices

If A is an n x m matrix, we can ask how many submatrices of order k are in A. Recall that we are only considering submatrices formed by consecutive rows and columns. There are (n – k + 1) sequences of consecutive rows of length k, such as 1:k, 2:(k+1), and so forth. Similarly, there are (m – k + 1) sequences of consecutive columns of length k. So the following SAS/IML function counts the number of submatrices of order k. You can call the function for different k=1, 2, 3, and 4 in order to solve the "Counting Squares" brain teaser:

/* count all submatrices of order k for the matrix A */ start CountAllSubmat(A, k); numRows = nrow(A)-k+1; numCols = ncol(A)-k+1; return numRows * numCols; finish; n = nrow(A); numSubmat = j(n, 1, .); do k = 1 to n; numSubmat[k] = CountAllSubmat(A, k); end; order = T(1:n); print order numSubmat; |

If you add up the number of squares each size, you get a total of 30 squares.

### Iterate over submatrices of a specified order

You can iterate over all submatrices of a specified order. You can perform a matrix computation on each submatrix (for example, compute its determinant), or you can pack it into a list for future processing. The following SAS/IML function iterates over submatrices and puts each submatrix into a list. The STRUCT subroutine from the ListUtil package can be used to indicate the contents of the list in a concise form.

start GetListSubmat(A, k); numSub = CountAllSubmat(A, k); L = ListCreate( numSub ); /* create list with numSub items */ cnt = 1; do i = 1 to nrow(A)-k+1; /* for each row */ do j = 1 to ncol(A)-k+1; /* and for each column */ B = submat(A, i, j, k); /* compute submatrix of order k from (i,j)th cell */ call ListSetItem(L, cnt, B); /* put in list (or compute with B here) */ cnt = cnt + 1; end; end; return L; finish; /* Example: get all matrices of order k=2 */ L = GetListSubmat(A, 2); package load ListUtil; call Struct(L); |

The output from the STRUCT call shwos that there are nine 2 x 2 matrices. The Value1–Valuie4 columns show the values (in row-major order). For example, the first 2 x 2 submatrix (in the upper left corner of A) is {4 3, 2 4}. The second submatrix (beginning with A[1,2]) is {3 1, 4 3}.

It is straightforward to modify the function to compute the determinant or some other matrix computation.

### Summary

This article shows some tips and techniques for dealing with submatrices of a matrix. The submatrices in this article are formed by using consecutive rows and columns. Thus, they are similar to the "Counting Squares" brain teaser, which requires that you consider sub-squares of order 1, 2, 3, and so forth.