Inside CNNs: Understanding Computer Vision Through Matrix Operations

Written by Aksel Aghajanyan · Aqwel AI Research.

---

Abstract

Computer vision systems do not “see” photographs the way humans do—they accumulate evidence through arithmetic on grids. At the machine level, an image is a structured grid of numbers—typically arranged as matrices or higher-order tensors—and nearly every practical pipeline, from classical filters to convolutional neural networks (CNNs), manipulates those numbers through linear algebra and localized reductions. This article explains how grayscale and RGB images become matrices, how discrete convolution slides kernels across spatial support to produce feature maps, how pooling downsamples activations, and why these operations map efficiently to modern GPUs. We walk through concrete numerical examples, NumPy implementations, and engineering considerations that practitioners encounter when building or debugging vision models. The goal is not to catalogue every architecture variant, but to build durable intuition for what tensors represent, how dimensions evolve layer by layer, and why matrix-centric computation remains the backbone of visual AI.

We deliberately emphasize reproducible micro-examples—small matrices where every arithmetic step can be checked by hand or in a notebook—because that discipline transfers directly to diagnosing shape mismatches and numeric anomalies in production stacks.

---

1. Introduction

Modern convolutional networks owe their dominance in perception tasks to a simple structural idea: reuse the same small set of weights across spatial positions. That reuse turns naturally into repeated inner-product-like operations between patches of an input tensor and filter weights—operations that hardware can parallelize aggressively. Before discussing depth, residual shortcuts, or attention mechanisms, it is worth grounding the reader in the discrete objects those networks consume and emit: tensors whose axes encode batch, channels, height, and width.

This document is written from an engineering lab perspective. When we say “convolution” in deep learning, we usually mean cross-correlation with a flipped kernel or, more commonly, cross-correlation without flipping—frameworks differ, but the essential mechanics—sliding window, multiply-accumulate, stride, padding—remain the same for intuition. We will be explicit about shapes, indexing, and boundary effects because those details determine memory use, receptive fields, and numerical stability in production systems.

Readers should leave with operational answers to debugging questions such as: “Why did my activation tensor shrink unexpectedly?” (check kernel size, stride, padding), “Why is VRAM exploding after I doubled resolution?” (activations scale roughly with pixel count per layer), and “Why does exporting ONNX change layouts?” (memory format versus logical axes). Those questions arise weekly in applied teams shipping vision models.

---

Reading Map

---

2. What Is an Image in Computer Vision?

A digital image is a sampling of a continuous scene onto a finite grid. Each sample location is a pixel (picture element). In software, pixels are stored as integers or floating-point numbers according to a convention:

The phrase “image file” (PNG, JPEG) describes encoding; once decoded into an array, the image is just numbers with shape (height, width) or (height, width, channels).

---

3. Images as Matrices

For grayscale, we can view an image as a matrix I with H rows and W columns:

I ∈ ℝ^(H×W)

Entry I[r, c] stores the intensity at row r, column c. Row-major versus column-major layout is a storage detail (NumPy defaults to row-major C-style); what matters for algorithms is consistent indexing.

For mini-batches in deep learning, we add axes:

• NCHW (batch, channels, height, width)—common in PyTorch defaults for conv modules.
• NHWC—often used in TensorFlow historically.

The choice affects stride patterns in memory but not the mathematical definition of convolution.

---

4. Grayscale Matrix Visualization

Consider a toy 3×3 grayscale patch:

[
 [  0,  50, 120 ],
 [ 80, 200, 255 ],
 [ 30, 100, 180 ]
]

Interpretation:

• Darker regions correspond to smaller values (closer to 0).
• Brighter regions correspond to larger values (closer to 255 in 8-bit displays).
• Spatial structure is preserved: neighbors on the grid remain neighbors in the scene (modulo lens distortion and projection).

When rendered to a screen, mapping [0,255] through a gamma curve produces perceptually uniform brightness for humans; CNNs instead learn filters directly on stored intensities or normalized tensors.

---

5. Why Black ≈ 0 and White ≈ 255 (for 8-bit)

Eight-bit channels afford 256 distinct levels—enough for consumer displays and compact storage. Setting black to 0 and white to 255 is a convention, not physics. Scientific imaging may use 16-bit depth or floating HDR; ML pipelines often rescale to [0,1] or standardize per dataset ((x − μ) / σ). The engineering takeaway: always know your preprocessing, because filter responses scale linearly with input magnitude unless normalized.

---

6. RGB Image Representation

Color images stack three channels—Red, Green, Blue—per pixel. The natural matrix generalization is a tensor X with shape (3, H, W) in channel-first layout:

X[k, r, c] = intensity of channel k at spatial location (r, c)

Human perception does not weight channels equally; camera sensors differ; illumination shifts absolute triplets. CNNs learn channel-mixing weights in the first convolutional layers (and beyond), effectively discovering task-specific color transformations instead of relying on hand-crafted rules.

---

7. Matrix Operations in Computer Vision

Classical vision used explicit linear transforms:

• Shifting / affine transforms compose as matrix multiplications on homogeneous coordinates.
• Gaussian smoothing is a convolution with a normalized kernel.
• Edge detection convolves with finite-difference operators (e.g., Sobel).

CNNs generalize this idea: instead of fixing kernels, we learn them from data subject to architectural constraints (locality via small kernels, translation-equivariance when weights are shared).

Important distinctions:

---

8. What Is Convolution (in Deep Learning)?

Given a 2-D input patch X and kernel K (also called filter weights), each output element is a sum of products over spatial support. For un-flipped cross-correlation—what most frameworks implement—the window aligns kernel entries with input entries under it:

y[i,j] = Σ Σ X[i+u, j+v] · K[u,v]

With stride S, we advance the window by S pixels along each axis. With padding P, we embed the input inside a larger canvas (commonly zero-padding) so border pixels participate symmetrically.

8.1 Stride in Practice

Stride controls how far the kernel advances after each output location is computed. Larger stride reduces spatial resolution of the output without introducing extra learned parameters—important when designing networks that must run at video frame rates or on embedded accelerators. When stride equals kernel size with non-overlapping windows, the operation resembles non-overlapping patch extraction followed by linear transforms.

8.2 Padding Modes

Valid convolution computes outputs only where the kernel fully fits inside the input—output shrinks unless special handling is applied. Same padding (approximately) seeks to preserve spatial size when stride is 1 and kernel dimensions are odd by distributing pad evenly (implementation details vary slightly across frameworks). Reflect padding copies boundary pixels symmetrically—sometimes preferred to reduce edge artifacts in segmentation.

8.3 Output Size Formula (Training-Time Sanity Check)

For height axis (width is analogous):

H_out = floor((H_in + 2P − kH) / S) + 1

When batching multiple images, tensor ranks expand: (N, C_in, H, W) convolved with (C_out, C_in, kH, kW) kernels yields (N, C_out, H_out, W_out) before bias and activation.

8.4 Multi-Channel Convolution as Dot Products

Each output channel f combines all input channels:

Y[f, i, j] = Σ_c Σ_u Σ_v X[c, i+u, j+v] · W[f, c, u, v] + b[f]

Thus a convolution layer linearly mixes channels within each spatial neighborhood—early layers often implement approximate luminance–chrominance decorrelation learned from data rather than hand-coded transforms.

---

9. Convolution Kernel Visualization (Conceptual)

Imagine a 3×3 input and a 2×2 kernel:

Input matrix:

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

Kernel:

[
 [ 1,  0 ],
 [ 0, -1 ]
]

Stride 1, valid convolution (no padding): the top-left window covers:

[1, 2]
[4, 5]

Element-wise multiply with K:

1·1 + 2·0 + 4·0 + 5·(−1) = 1 − 5 = −4

Slide right one column:

[2, 3]
[5, 6]

Response:

2·1 + 3·0 + 5·0 + 6·(−1) = 2 − 6 = −4

Continuing yields a 2×2 output feature map:

[
 [ -4, -4 ],
 [ -4, -4 ]
]

Completing all four placements explicitly

Positions (i,j) for top-left corners of the kernel inside the 3×3 valid grid are (0,0), (0,1), (1,0), (1,1).

• Window (0,0):

[1, 2]
[4, 5]

Dot with K: 1 − 5 = −4.

• Window (0,1):

[2, 3]
[5, 6]

Dot: 2 − 6 = −4.

• Window (1,0):

[4, 5]
[7, 8]

Dot: 4 − 8 = −4.

• Window (1,1):

[5, 6]
[8, 9]

Dot: 5 − 9 = −4.

For this toy kernel and ascending ramp input, responses coincide; real imagery breaks symmetry—edges and textures yield heterogeneous maps. The enduring lesson is mechanical: multiply, add, slide.

Strided Example (Stride S = 2, Same Input and Kernel, Valid)

Evaluating only positions where the kernel still fits after advancing by 2 along columns skips overlapping placements; along rows similarly. For H = W = 3, k = 2, valid convolution with S = 2 typically yields a 1×1 output (depending on alignment conventions)—demonstrating why practitioners carefully pad high-resolution tensors before aggressive striding.

---

10. Edge Detection Example

Discrete derivatives approximate brightness gradients. A simple horizontal edge filter:

[
 [-1, 0, 1],
 [-2, 0, 2],
 [-1, 0, 1]
]

(Sobel-like structure.) Large absolute responses align with rapid intensity transitions—edges that downstream layers can compose into textures and parts.

---

11. Feature Maps Explained

Each convolutional layer produces a stack of feature maps—one per output channel. Each map answers: “Where does this learned pattern activate?” Early layers often resemble edge and blob detectors; deeper layers encode progressively abstract cues conditional on dataset and objective.

Shape dynamics (single input channel, simplified):

Input:    (H, W)
Kernel:   (kH, kW)
Stride S, padding P

Output height:
  H_out = floor((H + 2P − kH) / S) + 1

Output width:
  W_out = floor((W + 2P − kW) / S) + 1

Wheninput has C_in channels and the layer has C_out filters, weights form a tensor of shape (C_out, C_in, kH, kW) in PyTorch Conv2d.

---

12. Pooling Operations

Pooling reduces spatial resolution while retaining dominant signals within neighborhoods.

A 2×2 max pool with stride 2 halves height and width (subject to divisibility). Pooling introduces translation tolerance locally and expands effective receptive fields of deeper layers without adding parameters at that operation.

12.1 Pooling as Downsample Without Trainable Weights

Unlike convolution, classical pooling layers introduce no new learnable scalar weights (aside from architectural variants). That makes them cheap at inference but also rigid—modern networks sometimes replace pooling with strided convolution to learn downsampling jointly with feature extraction.

12.2 NumPy: Max Pool 2×2, Stride 2

import numpy as np

def max_pool2x2(X: np.ndarray) -> np.ndarray:
    H, W = X.shape
    assert H % 2 == 0 and W % 2 == 0
    Y = np.zeros((H // 2, W // 2), dtype=X.dtype)
    for i in range(0, H, 2):
        for j in range(0, W, 2):
            Y[i // 2, j // 2] = np.max(X[i : i + 2, j : j + 2])
    return Y


pooled = max_pool2x2(np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]))
print(pooled)

---

13. How CNN Layers Work Together

A minimal CNN stacks:

1. Convolution — localized filtering + channel mixing.
2. Nonlinearity (e.g., ReLU) — introduces compositionality; without it, stacked linear ops collapse to a single linear map per spatial position up to degeneracies.
3. Pooling or strided convolution — downsampling.
4. Repeat.
5. Global pooling or flatten + linear — map to task logits.

The receptive field grows with depth and pooling—even neurons deep in the network indirectly depend on increasingly large input regions.

---

14. Why Matrices (and Tensors) Are Powerful in AI

• Structured locality: Nearby pixels correlate; convolution exploits that prior cheaply.
• Parameter efficiency: Sharing kernels across positions avoids separate weights per coordinate.
• Hardware alignment: Multiply-accumulate arrays on GPUs and TPUs match convolution inner loops.
• Composability: Depth stacks nonlinear feature hierarchies.

---

15. Engineering Perspective: GPUs and Tensor Cores

GPUs achieve throughput by scheduling many independent MAC lanes. Convolutions decompose into batched matrix multiplications via im2col or specialized fused kernels; tensor cores accelerate reduced-precision GEMMs used in mixed-precision training. From a systems viewpoint, memory bandwidth often bounds convolution-heavy workloads—another reason small kernels and efficient layouts matter.

---

16. Real-World Applications

16.1 Engineering Notes per Domain

Automotive: Camera clocks, rolling shutter, and HDR pipelines alter raw tensors before any neural net sees them—calibration matrices rectify projections; temporal fusion stacks tensors across time.

Medical imaging: Hounsfield units (CT) or MR intensities use different scales than RGB—normalization layers and modality-specific augmentations dominate reproducibility.

OCR: Text appears at arbitrary scales; multi-scale CNN features or transformer decoders align sequences of predictions with character or token alphabets.

Robotics: Depth sensors produce aligned RGB-D tensors where convolution may operate separately per modality before fusion layers.

---

17. Limitations and Challenges

• Noise sensitivity: Imaging sensors and compression artifacts perturb pixel values.
• Adversarial examples: Small deliberate perturbations can flip predictions despite imperceptible visual change.
• Compute & memory: High-resolution inputs explode activation tensor sizes.
• Data dependence: Biased datasets yield biased detectors.
• Annotation cost: Supervised training demands labeled examples at scale.

17.1 Mitigations Teams Actually Ship

• Input diversification: heavy augmentations (photometric, geometric) raise robustness—implemented as fast GPU-side tensor ops.
• Adversarial training: inject worst-case bounded noise during optimization—expensive but effective for security-sensitive deployments.
• Efficient architectures: depthwise separable convolutions, bottleneck blocks, and neural architecture search trade accuracy for FLOPs.
• Dataset audits: measure performance across demographic or environmental slices before release.

---

18. Future Directions (Brief)

Hybrid architectures combine convolutions with attention; Vision Transformers patchify images and apply self-attention across sequences of embeddings. Engineering trade-offs remain: inductive biases of CNNs versus flexibility of attention; latency versus accuracy; on-device quantization versus cloud-scale training.

18.1 From Convolution to Attention (Orientation Only)

Convolution assumes stationary local statistics—reasonable for textures and parts that recur across an image. Self-attention allows global interactions between all patch embeddings at a layer, trading compute (O(n²) in patches) for flexibility. Vision Transformers (ViT) reshape an image into a sequence of flattened patches, linearly embed each patch, add positional encodings, then stack transformer encoder blocks. CNNs remain dominant when data are limited or latency budgets are tight; ViT-class models shine at scale with ample supervision or strong pretraining.

18.2 Residual Connections (Conceptual)

Deep plain CNNs suffer optimization difficulties as depth grows. Residual networks introduce identity shortcuts y = F(x) + x, enabling gradients to bypass layers when F learns near-zero mappings—effectively letting depth be available without forcing every layer to distort representations catastrophically early in training.

---

Supplemental Technical Notes

S1. Why Convolution Is Efficient Compared With Fully Connected Layers on Images

A fully connected layer mapping H·W·C_in inputs to H·W·C_out outputs naïvely requires enormous weight matrices and ignores spatial locality. Convolution exploits:

1. Local connectivity: each output depends only on a k×k neighborhood per channel group.
2. Weight sharing: the same kernel applies at every spatial coordinate.

Parameter count scales with kernel area and channel depths—not directly with image width and height—making high-resolution inputs tractable.

S2. im2col and GEMM

Libraries often transform convolution into matrix multiplication by unfolding input patches into columns (im2col) and multiplying by a reshaped kernel matrix. That aligns convolution inner loops with highly tuned GEMM kernels on GPUs—another reason vision workloads benefit from tensor cores optimized for mixed-precision GEMMs.

S3. Tensors Beyond Four Dimensions

Video introduces a time axis; volumetric medical scans use depth spatial axis; batches stack N examples. The mental model remains identical: identify which axes participate in localized mixing (convolution across chosen spatial axes) versus global mixing (attention or FC).

S4. Edge Cases Practitioners Hit

• Even-sized kernels with “same” padding require asymmetric pad splits (floor/ceil)—verify framework behavior when porting weights between libraries.
• Dilated convolution expands receptive fields without pooling by inserting gaps in kernel sampling grids—common in segmentation decoders.
• Grouped convolution partitions channels to reduce parameter count (used in depthwise separable designs).

S5. Inner Products, Matrices, and Why GPUs Win

Each convolution response at a fixed location is an inner product between a flattened patch vector and a flattened kernel vector (plus bias). Stacking many patch vectors as rows of a matrix allows one large multiply—this is the algebraic heart of im2col. General matrix–matrix multiply (GEMM) saturates wide SIMD units and tensor cores because operands are dense and regular—convolution inherits that benefit after rewriting.

S6. Numerical Precision in Production

Training often uses float32 or mixed float16 with loss scaling; inference on edge devices may quantize to int8. Quantized convolution replaces real multiplies with integer equivalents plus scaling factors—smallkernels amplify quantization error at borders; calibration datasets tune per-channel scales.

When exporting models for mobile NPUs, verify that padding modes, rounding rules for integer divisions in pooling, and fused activation ordering match training graphs—silent divergence there yields tensors that look plausible yet degrade accuracy by whole percentage points.

---

19. Conclusion

Computer vision at the implementation layer is tensor arithmetic with clear rules: indexing, padding, stride, nonlinearity, pooling. Understanding discrete convolution demystifies feature maps and explains why CNNs parallelize effectively. For practitioners, the payoff is practical: informed choices about input resolution, kernel sizes, depth, and normalization prevent costly blind spots during training and deployment.

Long-term progress will blend classical convolutional priors with flexible sequence models where datasets and budgets justify them; nevertheless, matrix/tensor manipulations remain the substrate—whether expressed as explicit loops in NumPy, fused kernels in CUDA, or autograd graphs in PyTorch. Maintaining literacy in those foundations separates robust engineering from brittle integration.

---

References

1. LeCun, Y., Bengio, Y., & Hinton, G. (2015). Deep learning. Nature, 521(7553), 436–444.
2. Dumoulin, V., & Visin, F. (2016). A guide to convolution arithmetic for deep learning. arXiv:1603.07285.
3. Paszke, A., et al. (2019). PyTorch: An imperative style, high-performance deep learning library. NeurIPS.
4. Olah, C., et al. Distill.pub articles on convolutions and feature visualization (online).

---

Recommended Further Reading

• CS231n lecture notes (Stanford): Convolutional Neural Networks chapter.
• Deep Learning (Goodfellow, Bengio, Courville): convolution and structured probabilistic models sections.
• NVIDIA CUDA Deep Learning Developer Guide: sections on tensor layouts and convolution algorithms (implicit GEMM).

---

GitHub Project Ideas

1. Implement conv2d_naive_numpy and benchmark against scipy.signal.convolve2d modes.
2. Visualize receptive fields with gradient-based attribution on a small CNN.
3. Train a shallow CNN on MNIST entirely in NumPy + manual backprop for education.
4. Build a minimal im2col function for a single-channel image and compare outputs against a reference conv2d (shape checks only—speed is not the goal yet).

---

Future Research Directions (Aqwel AI Lab Notes)

• Differentiable visualization bridges between learned filters and symbolic explanations.
• Efficient convolution variants for scientific imaging modalities beyond RGB.
• Robustness certificates against bounded pixel perturbations in safety-critical sensing.
• Benchmark suites that pair tensor-level assertions (shape, dtype, stride) with pixel-level invariants for regression testing during refactors.

---

Appendix A — NumPy: Matrix Creation and Display

import numpy as np

# Toy grayscale patch (intensity 0–255)
I = np.array(
    [
        [0, 50, 120],
        [80, 200, 255],
        [30, 100, 180],
    ],
    dtype=np.float32,
)

print("Patch shape:", I.shape)
print(I)

---

Appendix B — Matplotlib: Visualize as Image

import matplotlib.pyplot as plt

# Normalize for display if needed
plt.imshow(I, cmap="gray", vmin=0, vmax=255)
plt.colorbar(label="Intensity")
plt.title("Grayscale patch as an image")
plt.axis("off")
plt.show()

---

Appendix C — Convolution from Scratch (Valid, stride 1)

import numpy as np


def conv2d_valid(X: np.ndarray, K: np.ndarray) -> np.ndarray:
    """2-D valid convolution (cross-correlation) with stride 1, no padding."""
    H, W = X.shape
    kH, kW = K.shape
    out_h = H - kH + 1
    out_w = W - kW + 1
    Y = np.zeros((out_h, out_w), dtype=X.dtype)
    for i in range(out_h):
        for j in range(out_w):
            region = X[i : i + kH, j : j + kW]
            Y[i, j] = np.sum(region * K)
    return Y


X = np.arange(1, 10, dtype=np.float32).reshape(3, 3)
K = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.float32)
print("Input:\n", X)
print("Kernel:\n", K)
print("Output:\n", conv2d_valid(X, K))

---

Appendix D — Edge Emphasis with OpenCV (Optional)

# pip install opencv-python  (optional dependency)
import cv2
import numpy as np

img = np.uint8(I)  # from Appendix A
gx = cv2.Sobel(img, cv2.CV_32F, 1, 0, ksize=3)
gy = cv2.Sobel(img, cv2.CV_32F, 0, 1, ksize=3)
mag = np.sqrt(gx * gx + gy * gy)

---

Appendix E — Global Average Pooling Concept

Replace large spatial maps with channel-wise averages:

For channel c:  GAP_c = (1 / (H·W)) Σ_{r,c} A[c, r, c]

This reduces parameters compared with large fully connected stacks and helps localization tasks.

---

Appendix F — Matplotlib: Synthetic Feature Map Heatmap

import matplotlib.pyplot as plt
import numpy as np

rng = np.random.default_rng(0)
H, W = 8, 10
# toy activation map after conv — structured + noise
r = np.linspace(-1, 1, H)
c = np.linspace(-1, 1, W)
R, C = np.meshgrid(c, r)
feature_map = np.exp(-(R ** 2 + C ** 2)) + 0.15 * rng.standard_normal((H, W))

plt.figure(figsize=(7, 5))
plt.imshow(feature_map, cmap="magma", interpolation="nearest")
plt.colorbar(label="Activation")
plt.title("Toy convolution response (illustrative)")
plt.xlabel("Width index")
plt.ylabel("Height index")
plt.tight_layout()
plt.show()

---

Appendix G — Diagnostics Practitioners Run

• Activation histograms: detect dying ReLUs or exploding norms layer-wise.
• Gradient norms: identify unstable depths before adding normalization or residual wiring.
• Kernel visualizations: first-layer filters often approximate Gabor-like edge detectors when trained on natural images—sanity check that preprocessing is correct.

---

Visual Assets to Commission (Design Brief)

---

Optional Advanced Note — PyTorch Tensors

PyTorch represents batches as torch.Tensor objects with dtype and device (cpu / cuda). Autograd tracks operations for reverse-mode differentiation—what makes end-to-end CNN training feasible at scale.

---

Technical note: This article uses discrete summation notation appropriate for implemented CNNs. Continuous-time convolution integrals from signal processing provide historical context but are not required to implement modern vision networks.