Covariance Matrix from Image Patches

hardFREEVectorsCV: Mathematical Foundations

Compute the covariance matrix of image patches for PCA-based texture analysis.

Given an image and patch size, extract all overlapping patches, flatten them to vectors, and compute the covariance matrix. This is a fundamental step in texture analysis and eigenface computation.

Algorithm:

Extract all overlapping patches of size $p \times p$ from the image
Flatten each patch (row by row) into a vector of length $p^2$
Stack vectors as rows of data matrix $X$ (shape: n_patches × p²)
Center the data: $X_c = X - \bar{X}$ (subtract mean vector)
Compute covariance: $\Sigma = \frac{1}{n-1} X_c^T X_c$

If there is only one patch ( $n = 1$ ), return a $p^2 \times p^2$ matrix of zeros.

The covariance matrix $\Sigma$ captures correlations between pixel positions within patches.

Applications:

PCA for dimensionality reduction
Eigenface/Fisherface recognition
Texture synthesis and analysis

Example

Input

image = [[1, 2, 3],
         [4, 5, 6],
         [7, 8, 9]]
patch_size = 2

Output

[[3.3333, 3.3333, 3.3333, 3.3333],
 [3.3333, 3.3333, 3.3333, 3.3333],
 [3.3333, 3.3333, 3.3333, 3.3333],
 [3.3333, 3.3333, 3.3333, 3.3333]]

Reasoning

Input: a $3 \times 3$ grayscale image and $p = 2$ . A $2 \times 2$ window fits in $3 - 2 + 1 = 2$ positions down and 2 across, so there are $n = 4$ overlapping patches.
Patches as vectors: reading each window row by row gives $(1,2,4,5)$ , $(2,3,5,6)$ , $(4,5,7,8)$ and $(5,6,8,9)$ . These are the 4 rows of $X$ , which is $4 \times 4$ because $p^2 = 4$ .
Centering: the mean of each column (pixel position) is $(3, 4, 6, 7)$ . Subtracting it, every column of $X_c$ becomes $(-2, -1, 1, 2)^T$ .
Covariance: entry $\Sigma_{ij}$ is the dot product of centred columns $i$ and $j$ , divided by $n - 1 = 3$ : $\frac{(-2)^2 + (-1)^2 + 1^2 + 2^2}{3} = \frac{10}{3} \approx 3.3333$ for every $i, j$ .
Output: a $p^2 \times p^2 = 4 \times 4$ symmetric matrix, rounded to 4 decimal places. Notice it is constant: in this linear ramp all four pixel positions rise and fall together, so every variance equals every covariance.

Constraints

image: 2D numpy array (grayscale)
patch_size: integer p (patches are p×p)
Return: covariance matrix of shape (p², p²)
Round to 4 decimal places

Learn this in

Mathematical FoundationsFoundations

More from CV: Mathematical Foundations

Adam Optimizer Stephard Power Iteration for Dominant Eigenvectorhard SVD Image Compressionhard Batch Affine Transform with PyTorchmedium Gram-Schmidt Orthonormalizationmedium

Concept primer

Background Knowledge

Covariance matrices are fundamental statistical tools that capture the linear relationships between variables in a dataset. In the context of image patches, each pixel position within a patch is treated as a variable, and the covariance matrix reveals which pixel positions tend to vary together. When you flatten a patch into a vector, you're converting spatial information into a feature space where the covariance matrix becomes a symmetric, positive semi-definite matrix of shape $p^2 \times p^2$ .

Principal Component Analysis (PCA) relies on diagonalizing the covariance matrix to identify the directions of maximum variance in your data. For image patches, this means finding the most important patterns or textures that explain variation across your dataset. The eigenvectors of the covariance matrix form an orthonormal basis, and the eigenvalues tell you how much variance is explained along each direction. This is why computing the covariance matrix accurately is critical—it directly determines which features PCA will extract.

Centering the data is essential because the covariance matrix measures deviations from the mean. Without centering, you're computing a second moment rather than a true covariance, which would incorrectly include information about the mean intensity of patches rather than their internal structure. The factor of $\frac{1}{n-1}$ (Bessel's correction) provides an unbiased estimate when working with sample data rather than an entire population.

Algorithm/Approach

The solution follows a data preprocessing pipeline common in machine learning:

Spatial extraction: Slide a $p \times p$ window across the image with stride 1 to capture all overlapping patches
Vectorization: Convert 2D patches to 1D vectors to create a data matrix where each row is a flattened patch
Normalization: Center the data by subtracting the mean vector computed across all patches
Covariance computation: Apply the mathematical formula $\Sigma =$ \frac{1}{n-1} $X_c^T X_c$ to obtain the final matrix

This approach is widely used in texture analysis and face recognition systems, where local patch statistics reveal meaningful patterns.

Step-by-Step Strategy

Step 1: Extract Overlapping Patches

Iterate through the image with a sliding window of size $p \times p$
For an image of size $H \times W$ , you'll have $(H - p + 1) \times (W - p + 1)$ patches
Store each patch as a separate entity (you can use a list or pre-allocate an array)

Step 2: Flatten and Stack into Data Matrix

Convert each $p \times p$ patch into a 1D vector of length $p^2$
Stack all vectors as rows in a matrix $X$ of shape $n$ \text{ patches} $\times p^2$
Ensure you're working with floating-point data types to avoid precision loss

Step 3: Compute the Mean Vector

Calculate $\bar{X}$ as the mean of each column (across all patches)
This represents the "average pixel value" at each position within patches

Step 4: Center the Data

Subtract the mean vector from every row: $X_c = X - \bar{X}$
This ensures the covariance captures variance, not absolute intensity

Step 5: Compute Covariance Matrix

Use matrix multiplication: $\Sigma =$ \frac{1}{n-1} $X_c^T X_c$
The result is a $p^2 \times p^2$ symmetric matrix
Verify symmetry as a sanity check (numerical errors may cause small asymmetries)

Common Pitfalls

Forgetting to center: Computing $X^T X$ without centering gives incorrect results that conflate mean and variance
Off-by-one errors in patch extraction: Ensure your loop bounds correctly capture all overlapping patches without duplicates or gaps
Data type issues: Using integer types for pixel values can cause overflow or precision loss; convert to float early
Memory efficiency: For large images and small patch sizes, the number of patches can be enormous. Consider whether you need to store all patches or can compute statistics incrementally
Normalization factor confusion: Remember it's $\frac{1}{n-1}$ for sample covariance (unbiased), not $\frac{1}{n}$ (which is biased)
Matrix dimension mismatch: Verify that $X_c$ has shape $n \times p^2$ before computing $X_c^T X_c$ , which should yield $p^2 \times p^2$

Time & Space Complexity

Time Complexity:

Patch extraction: $O((H - p + 1)(W - p + 1) \cdot p^2)$ to extract and flatten all patches
Centering: $O(n \cdot p^2)$ where $n = (H - p + 1)(W - p + 1)$
Covariance computation: $O(n \cdot p^4)$ for the matrix multiplication $X_c^T X_c$
Overall: $O(n \cdot p^4)$ where the dominant cost is the matrix multiplication

Space Complexity:

Data matrix $X$ : $O(n \cdot p^2)$
Covariance matrix $\Sigma$ : $O(p^4)$
Overall: $O(n \cdot p^2 + p^4)$

For typical values (e.g., $p = 8$ , image $256 \times 256$ ), this is manageable, but for larger patches or images, memory becomes a constraint. Some implementations use incremental covariance computation to reduce memory usage.

solution.py

Python 3.13

Test results