Matrix Operations

Beginner: matrix addition is straightforward — add matching cells, and it only works when both matrices are exactly the same shape. Matrix multiplication is the one to understand deeply: to get one output number, take a full row of the first matrix, a full column of the second, and compute their dot product.

Intermediate: the shape rule is worth memorizing precisely, because it's the single most common source of bugs when building neural networks by hand: multiplying an (m×n) matrix by an (n×p) matrix requires the inner dimensions to match (both n), and the result is an (m×p) matrix — the outer dimensions "survive." If the inner dimensions don't match, the multiplication is simply undefined, not approximately correct.

Advanced: matrix multiplication composes transformations. If matrix B rotates space and matrix A then stretches it, the single matrix AB does both at once, in that order, to any vector you feed it. This is exactly how a multi-layer neural network's effective transformation could in principle be written as one giant matrix — if there were no non-linear activation functions breaking up the chain (which is precisely why those non-linearities are essential: without them, an arbitrarily deep network would collapse to a single linear transformation).

This four-line snippet is a full forward pass through one neural network layer — W @ x + b is matrix multiplication plus a bias vector, exactly as described above, and np.maximum(0, z) is the ReLU non-linearity that keeps stacked layers from collapsing into a single linear transformation.

#include <algorithm>
#include <cstdio>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

std::vector<double> matVecMul(const Matrix& M, const std::vector<double>& v) {
    std::vector<double> out(M.size(), 0.0);
    for (size_t i = 0; i < M.size(); ++i)
        for (size_t j = 0; j < v.size(); ++j)
            out[i] += M[i][j] * v[j];
    return out;
}

std::vector<double> addVec(const std::vector<double>& u, const std::vector<double>& v) {
    std::vector<double> out(u.size());
    for (size_t i = 0; i < u.size(); ++i) out[i] = u[i] + v[i];
    return out;
}

std::vector<double> relu(const std::vector<double>& v) {
    std::vector<double> out(v.size());
    for (size_t i = 0; i < v.size(); ++i) out[i] = std::max(0.0, v[i]);
    return out;
}

int main() {
    Matrix W = {{0.2, -0.5}, {0.8, 0.1}};   // weights, shape (2, 2)
    std::vector<double> x = {1.0, 2.0};       // input, shape (2,)
    std::vector<double> b = {0.1, -0.2};      // bias, shape (2,)

    std::vector<double> z = addVec(matVecMul(W, x), b);  // matmul then add bias
    std::vector<double> output = relu(z);                  // ReLU activation

    std::printf("[%.4f, %.4f]\n", output[0], output[1]);
    return 0;
}