Linear Algebra I: Vectors & Matrices

A vector has two equivalent faces, and you need both.

The geometric face: A vector is an arrow in space. It has a direction and a magnitude (length), but no fixed starting point — two arrows are the same vector if they point the same way and have the same length. In 2D this is an arrow you could draw on graph paper; in 3D it points somewhere in the room. This is the physicist's picture.

The algebraic face: A vector is an ordered tuple of numbers, written as a column: $$\mathbf{v} = \begin{pmatrix} v_1 \ v_2 \ \vdots \ v_n \end{pmatrix}$$ Each number $v_i$ is a component of v along the $i$-th coordinate axis.

Why both? The geometric view gives you intuition about what operations mean. The algebraic view gives you a way to compute with them. Throughout linear algebra you will constantly translate between these two faces.

Vector addition is done component-wise: $(\mathbf{u} + \mathbf{v})_i = u_i + v_i$. Geometrically, place the tail of v at the tip of u — the sum is the arrow from the tail of u to the tip of v. This is the parallelogram law, and it captures how forces combine in physics.

Scalar multiplication stretches (or shrinks, or flips) the arrow without changing its direction (unless the scalar is negative): $(c\mathbf{v})_i = c \cdot v_i$.

The norm (length) of a vector in $\mathbb{R}^n$ is the Euclidean distance from origin to tip: $$|\mathbf{v}| = \sqrt{v_1^2 + v_2^2 + \cdots + v_n^2}$$

A unit vector has norm 1. To convert any nonzero vector to a unit vector, divide by its norm: $\hat{\mathbf{v}} = \mathbf{v}/|\mathbf{v}|$.

import numpy as np

u = np.array([1.0, 2.0, 3.0])
v = np.array([4.0, -1.0, 0.0])

print(u + v)        # [5.  1.  3.]
print(3 * u)        # [3.  6.  9.]
print(np.linalg.norm(u))   # sqrt(14) ≈ 3.742
print(u / np.linalg.norm(u))  # unit vector of u

The dot product (also called inner product or scalar product) takes two vectors and returns a number: $$\mathbf{u} \cdot \mathbf{v} = \sum_{i=1}^{n} u_i v_i = u_1 v_1 + u_2 v_2 + \cdots + u_n v_n$$

But this algebraic formula hides a beautiful geometric meaning. The dot product equals: $$\mathbf{u} \cdot \mathbf{v} = |\mathbf{u}| |\mathbf{v}| \cos\theta$$ where $\theta$ is the angle between the two vectors.

What does $\cos\theta$ tell us?

• If $\theta = 0°$ (same direction): $\mathbf{u} \cdot \mathbf{v} = |\mathbf{u}||\mathbf{v}|$ (maximum positive)
• If $\theta = 90°$ (perpendicular): $\mathbf{u} \cdot \mathbf{v} = 0$ (orthogonal — one of the most important special cases)
• If $\theta = 180°$ (opposite directions): $\mathbf{u} \cdot \mathbf{v} = -|\mathbf{u}||\mathbf{v}|$ (maximum negative)

The dot product is literally a measure of alignment between two vectors. This is exactly why it appears everywhere in machine learning: the similarity between a query and a document in search, the activation in a neuron (weights dot input), the projection of one vector onto another direction.

Projection: The scalar projection of u onto v is $\text{proj}_{\mathbf{v}} \mathbf{u} = (\mathbf{u} \cdot \hat{\mathbf{v}})\hat{\mathbf{v}}$, which is the "shadow" of u cast onto the line through v.

Cauchy-Schwarz inequality: $|\mathbf{u} \cdot \mathbf{v}| \leq |\mathbf{u}||\mathbf{v}|$. This is equivalent to $|\cos\theta| \leq 1$, but the algebraic proof works in any dimension and is one of the most important inequalities in all of mathematics.

import numpy as np

u = np.array([1.0, 0.0])
v = np.array([1.0, 1.0]) / np.sqrt(2)  # unit vector at 45 degrees

dot = np.dot(u, v)
angle_rad = np.arccos(dot / (np.linalg.norm(u) * np.linalg.norm(v)))
print(f"dot product: {dot:.4f}")          # 0.7071 = cos(45°)
print(f"angle (degrees): {np.degrees(angle_rad):.1f}")  # 45.0

# Check orthogonality
a = np.array([1.0, 2.0, -1.0])
b = np.array([3.0, 0.0, 3.0])
print("orthogonal:", np.isclose(np.dot(a, b), 0))  # True

So far we have worked in $\mathbb{R}^n$, but the real power of linear algebra comes from the realization that many mathematical objects share the same structural skeleton.

Definition. A vector space over $\mathbb{R}$ is a set $V$ together with two operations — addition ($V \times V \to V$) and scalar multiplication ($\mathbb{R} \times V \to V$) — satisfying eight axioms:

1. Commutativity of addition: $\mathbf{u} + \mathbf{v} = \mathbf{v} + \mathbf{u}$
2. Associativity of addition: $(\mathbf{u}+\mathbf{v})+\mathbf{w} = \mathbf{u}+(\mathbf{v}+\mathbf{w})$
3. Zero vector: there exists $\mathbf{0}$ such that $\mathbf{v} + \mathbf{0} = \mathbf{v}$
4. Additive inverse: for every $\mathbf{v}$, there exists $-\mathbf{v}$ such that $\mathbf{v} + (-\mathbf{v}) = \mathbf{0}$
5. Scalar identity: $1 \cdot \mathbf{v} = \mathbf{v}$
6. Scalar associativity: $(ab)\mathbf{v} = a(b\mathbf{v})$ 7–8. Two distributive laws: $a(\mathbf{u}+\mathbf{v}) = a\mathbf{u} + a\mathbf{v}$ and $(a+b)\mathbf{v} = a\mathbf{v} + b\mathbf{v}$

Examples of vector spaces:

• $\mathbb{R}^n$ (the familiar one)
• The set of all polynomials of degree $\leq 3$: ${a_0 + a_1 x + a_2 x^2 + a_3 x^3}$
• The set of all continuous functions on $[0,1]$
• The set of all $m \times n$ real matrices

Non-examples: The positive reals under standard multiplication and addition do not form a vector space (no additive inverse). The set ${(x,y) : x^2 + y^2 = 1}$ (a circle) is not a vector space (not closed under addition).

Subspace: A subset $W \subseteq V$ is a subspace if it is itself a vector space, which reduces to three conditions: (1) $\mathbf{0} \in W$, (2) closed under addition, (3) closed under scalar multiplication. Every line or plane through the origin in $\mathbb{R}^3$ is a subspace.

# Verify subspace closure in R^3: the plane x + y + z = 0
# A subspace must pass through the origin and be closed under + and scalar mult.

def on_plane(v):
    return abs(v[0] + v[1] + v[2]) < 1e-10

u = [1.0, -1.0, 0.0]  # 1 - 1 + 0 = 0 ✓
v = [2.0, 0.0, -2.0]  # 2 + 0 - 2 = 0 ✓

print(on_plane([u[i] + v[i] for i in range(3)]))  # True: closed under addition
print(on_plane([3*u[i] for i in range(3)]))        # True: closed under scalar mult
print(on_plane([0.0, 0.0, 0.0]))                   # True: contains zero

These three concepts form the heart of the module. Get them right and everything else falls into place.

Linear combination: $c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2 + \cdots + c_k \mathbf{v}_k$ is a linear combination of vectors $\mathbf{v}_1, \ldots, \mathbf{v}_k$. The set of all linear combinations is the span of those vectors.

Linear independence: Vectors $\mathbf{v}_1, \ldots, \mathbf{v}_k$ are linearly independent if the only solution to $c_1 \mathbf{v}_1 + \cdots + c_k \mathbf{v}_k = \mathbf{0}$ is $c_1 = c_2 = \cdots = c_k = 0$. If any other solution exists, some vector is redundant — it lies in the span of the others.

Analogy: Think of building IKEA furniture. Linearly independent vectors are like instructions that each add genuinely new information. A dependent set is like having an instruction that says "do the same as steps 1 and 3 combined" — redundant.

Basis: A basis for a vector space $V$ is a set of vectors that is (1) linearly independent and (2) spans $V$. A basis is the most efficient, non-redundant set that can describe every vector in the space.

Dimension: Every basis for a given vector space has the same number of vectors, and this number is the dimension of the space, $\dim(V)$.

• Standard basis of $\mathbb{R}^n$: the $n$ unit vectors $\mathbf{e}_1 = (1,0,\ldots,0)^T$, etc.
• The space of polynomials of degree $\leq 3$ has dimension 4, with basis ${1, x, x^2, x^3}$.

Key fact: In $\mathbb{R}^n$, any $n$ linearly independent vectors form a basis. More than $n$ vectors must be dependent; fewer than $n$ vectors cannot span.

Coordinates: Once you fix a basis ${\mathbf{b}_1, \ldots, \mathbf{b}_n}$, every vector $\mathbf{v}$ in the space has a unique representation $\mathbf{v} = c_1 \mathbf{b}_1 + \cdots + c_n \mathbf{b}_n$. The scalars $(c_1, \ldots, c_n)$ are the coordinates of v in that basis. Different bases give different coordinates — the same vector, different descriptions.

import numpy as np

# Check if 3 vectors in R^3 are linearly independent
# by checking if the matrix they form has nonzero determinant
v1 = np.array([1.0, 0.0, 0.0])
v2 = np.array([0.0, 1.0, 0.0])
v3 = np.array([1.0, 1.0, 0.0])  # v3 = v1 + v2 -- dependent!

M = np.column_stack([v1, v2, v3])
print(f"det = {np.linalg.det(M):.4f}")  # 0.0000 -- linearly dependent

v3_indep = np.array([0.0, 0.0, 1.0])
M2 = np.column_stack([v1, v2, v3_indep])
print(f"det = {np.linalg.det(M2):.4f}")  # 1.0000 -- linearly independent (basis!)

A matrix is a rectangular array of numbers arranged in rows and columns. An $m \times n$ matrix $A$ has $m$ rows and $n$ columns, with entry $A_{ij}$ in row $i$, column $j$.

$$A = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \ a_{21} & a_{22} & \cdots & a_{2n} \ \vdots & & \ddots & \vdots \ a_{m1} & a_{m2} & \cdots & a_{mn} \end{pmatrix}$$

Two ways to see a matrix:

1. As a collection of column vectors: $A = [\mathbf{a}_1 \mid \mathbf{a}_2 \mid \cdots \mid \mathbf{a}_n]$ where each $\mathbf{a}_j \in \mathbb{R}^m$.
2. As a collection of row vectors: the $i$-th row $\mathbf{r}_i^T$ gives the $i$-th constraint or output.

Matrix-vector multiplication $A\mathbf{x}$ (for $A$ an $m\times n$ matrix and $\mathbf{x} \in \mathbb{R}^n$) is defined as: $$A\mathbf{x} = x_1 \mathbf{a}_1 + x_2 \mathbf{a}_2 + \cdots + x_n \mathbf{a}_n$$ This is a linear combination of the columns of $A$, weighted by the entries of x. This perspective is fundamental: $A\mathbf{x}$ sweeps out the column space of $A$ as x varies.

Alternatively, the $i$-th entry of $A\mathbf{x}$ is the dot product of the $i$-th row of $A$ with x: $(A\mathbf{x})_i = \mathbf{r}_i \cdot \mathbf{x}$.

Special matrices:

• Identity matrix $I_n$: ones on the diagonal, zeros elsewhere. $I\mathbf{x} = \mathbf{x}$ for all x.
• Zero matrix: all zeros.
• Diagonal matrix: nonzero entries only on the main diagonal. Scales each coordinate independently.
• Transpose $A^T$: $(A^T){ij} = A{ji}$ — flip rows and columns.
• Symmetric matrix: $A^T = A$. Covariance matrices in statistics are always symmetric.

import numpy as np

A = np.array([[1, 2, 3],
              [4, 5, 6]])  # 2x3 matrix
x = np.array([1.0, 0.0, -1.0])

# Matrix-vector product: linear combination of columns
result = A @ x
print(result)  # A[:,0]*1 + A[:,1]*0 + A[:,2]*(-1) = [1-3, 4-6] = [-2, -2]

# Verify via column perspective
col_combo = 1*A[:,0] + 0*A[:,1] + (-1)*A[:,2]
print(col_combo)  # same: [-2, -2]

print(A.T)  # transpose: shape becomes (3, 2)

Matrix multiplication $C = AB$ (where $A$ is $m \times k$ and $B$ is $k \times n$) produces an $m \times n$ matrix with: $$C_{ij} = \sum_{p=1}^{k} A_{ip} B_{pj} = \text{(row }i\text{ of }A) \cdot \text{(column }j\text{ of }B)$$

This looks mechanical, but there are four equivalent ways to see what's happening, each revealing something different:

1. Row $\times$ Column (standard): $C_{ij} = \mathbf{r}_i(A) \cdot \mathbf{c}_j(B)$. Standard definition.

2. Column view: Each column of $C$ is $A$ times the corresponding column of $B$: $\mathbf{c}_j(C) = A\mathbf{c}_j(B)$. So $AB$ applies $A$ to every column of $B$.

3. Row view: Each row of $C$ is the corresponding row of $A$ times $B$.

4. Sum of outer products: $AB = \sum_{p=1}^{k} \mathbf{a}_p \mathbf{b}_p^T$ where $\mathbf{a}_p$ is column $p$ of $A$ and $\mathbf{b}_p^T$ is row $p$ of $B$. Each term is a rank-1 matrix; $AB$ is a sum of $k$ rank-1 pieces. This view is crucial in deep learning (attention, factorization).

Key properties:

• Associative: $(AB)C = A(BC)$
• Distributive: $A(B+C) = AB + AC$
• NOT commutative in general: $AB \neq BA$ (in fact, the dimensions may not even match for $BA$)
• $(AB)^T = B^T A^T$ (order reverses under transpose)

Why is the order rule "inner dimensions must match"? Because $A$ maps $\mathbb{R}^k \to \mathbb{R}^m$ and $B$ maps $\mathbb{R}^n \to \mathbb{R}^k$. To compose them (apply $B$ first, then $A$), the output space of $B$ must equal the input space of $A$.

Computational cost: Naive matrix multiplication of two $n \times n$ matrices costs $O(n^3)$ operations. This is why matrix multiplication is the dominant operation in neural network training, and why efficient BLAS routines and GPU acceleration matter enormously.

import numpy as np

A = np.array([[1, 2], [3, 4], [5, 6]])  # 3x2
B = np.array([[7, 8, 9], [10, 11, 12]])  # 2x3

C = A @ B  # 3x3 result
print(C)
# Row 0 of C: [1*7+2*10, 1*8+2*11, 1*9+2*12] = [27, 30, 33]

# Non-commutativity: AB != BA in general
P = np.array([[1,2],[3,4]])
Q = np.array([[0,1],[1,0]])
print(P @ Q)   # [[2,1],[4,3]]
print(Q @ P)   # [[3,4],[1,2]]  -- different!

# Outer product sum (4th view)
col0_A = A[:, 0:1]  # column vector (3,1)
row0_B = B[0:1, :]  # row vector (1,3)
print(col0_A @ row0_B)  # rank-1 matrix

This concept is arguably the most important in the entire module. A matrix is not just an array of numbers — it is a linear transformation in disguise.

Definition. A function $T: V \to W$ is a linear transformation (or linear map) if for all vectors $\mathbf{u}, \mathbf{v}$ and scalars $c$:

1. $T(\mathbf{u} + \mathbf{v}) = T(\mathbf{u}) + T(\mathbf{v})$ (additivity)
2. $T(c\mathbf{v}) = cT(\mathbf{v})$ (homogeneity)

Equivalently: $T(c\mathbf{u} + d\mathbf{v}) = cT(\mathbf{u}) + dT(\mathbf{v})$ for all scalars $c, d$.

The key theorem: Every linear transformation $T: \mathbb{R}^n \to \mathbb{R}^m$ can be represented as multiplication by an $m \times n$ matrix. The matrix is built by seeing where $T$ sends each standard basis vector: column $j$ of the matrix is $T(\mathbf{e}_j)$.

Why does this work? Because every $\mathbf{x} = x_1\mathbf{e}_1 + \cdots + x_n\mathbf{e}_n$, so: $$T(\mathbf{x}) = x_1 T(\mathbf{e}_1) + \cdots + x_n T(\mathbf{e}_n) = [T(\mathbf{e}_1) \mid \cdots \mid T(\mathbf{e}_n)]\mathbf{x}$$

Examples in 2D (the matrix that represents each):

• Rotation by angle $\theta$: $\begin{pmatrix}\cos\theta & -\sin\theta \ \sin\theta & \cos\theta\end{pmatrix}$
• Reflection across x-axis: $\begin{pmatrix}1 & 0 \ 0 & -1\end{pmatrix}$
• Projection onto x-axis: $\begin{pmatrix}1 & 0 \ 0 & 0\end{pmatrix}$
• Shear: $\begin{pmatrix}1 & 1 \ 0 & 1\end{pmatrix}$

Four fundamental subspaces of a matrix $A$ ($m \times n$):

• Column space $C(A) \subseteq \mathbb{R}^m$: all possible outputs $A\mathbf{x}$; span of columns.
• Null space $N(A) \subseteq \mathbb{R}^n$: all x such that $A\mathbf{x} = \mathbf{0}$; inputs that vanish.
• Row space $C(A^T) \subseteq \mathbb{R}^n$: span of rows.
• Left null space $N(A^T) \subseteq \mathbb{R}^m$.

Rank-Nullity Theorem: $\text{rank}(A) + \dim(N(A)) = n$. The number of independent columns plus the number of dimensions that collapse to zero must equal the total input dimension. This is one of the most useful facts in linear algebra — it immediately tells you the structure of the solution set to $A\mathbf{x} = \mathbf{b}$.

import numpy as np

theta = np.pi / 4  # 45 degrees
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])

# Rotate the standard basis vectors
e1 = np.array([1.0, 0.0])
e2 = np.array([0.0, 1.0])
print(R @ e1)  # [0.707, 0.707] -- points 45 degrees from x-axis
print(R @ e2)  # [-0.707, 0.707] -- 90 degrees further

# Null space example: projection onto x-axis
P = np.array([[1.0, 0.0], [0.0, 0.0]])
v = np.array([3.0, 5.0])
print(P @ v)   # [3.0, 0.0] -- y-component collapses to 0
# Null space of P: all vectors [0, y] -- entire y-axis

The original motivation for all of linear algebra is solving systems like: $$\begin{cases} 2x_1 + x_2 - x_3 = 8 \ -3x_1 - x_2 + 2x_3 = -11 \ -2x_1 + x_2 + 2x_3 = -3 \end{cases}$$

In matrix form, this is $A\mathbf{x} = \mathbf{b}$. The system asks: does b lie in the column space of $A$? If so, which combination of columns of $A$ produces it?

Three possible outcomes (and their geometric meaning):

1. Unique solution: The columns of $A$ are linearly independent and b is in their span. Geometrically in 2D: two lines that intersect at exactly one point.
2. Infinitely many solutions: The columns are dependent and b happens to be in their span anyway. Geometrically in 2D: two lines that are actually the same line.
3. No solution: b is NOT in the column space of $A$. Geometrically in 2D: two parallel but distinct lines.

The augmented matrix $[A \mid \mathbf{b}]$ packs the entire system into one object for convenient manipulation.

Conditions for solvability:

• $A\mathbf{x} = \mathbf{b}$ has a solution if and only if $\text{rank}(A) = \text{rank}([A \mid \mathbf{b}])$.
• If $\text{rank}(A) = n$ (full column rank), the solution is unique.
• If $\text{rank}(A) < n$, the solution set is a translation of the null space: ${\mathbf{x}_p + \mathbf{x}_h : \mathbf{x}_h \in N(A)}$ where $\mathbf{x}_p$ is any particular solution.

The null space interpretation of solution structure is profound: the "freedom" in the solution (how many free parameters you have) is exactly the dimension of the null space — a purely intrinsic property of $A$ that has nothing to do with b.

import numpy as np

A = np.array([[2.,  1., -1.],
              [-3.,-1.,  2.],
              [-2., 1.,  2.]])
b = np.array([8., -11., -3.])

# Solve Ax = b
x = np.linalg.solve(A, b)
print("Solution:", x)          # [2. 3. -1.]
print("Verify Ax = b:", np.allclose(A @ x, b))  # True

# Check rank
print("rank(A):", np.linalg.matrix_rank(A))  # 3 (full rank -> unique solution)

Gaussian elimination is an algorithm for transforming a system $A\mathbf{x} = \mathbf{b}$ into an equivalent triangular system that is easy to solve by back substitution.

The key insight: you can perform three elementary row operations on the augmented matrix $[A \mid \mathbf{b}]$ without changing the solution set:

1. Row swap: swap row $i$ and row $j$
2. Row scaling: multiply row $i$ by nonzero scalar $c$
3. Row addition: add a multiple of row $j$ to row $i$

The algorithm (forward elimination):

1. Find the first column with a nonzero entry — this is the pivot column. Choose a pivot (preferably the largest entry for numerical stability — partial pivoting).
2. Use row additions to zero out all entries below the pivot.
3. Move to the next column/row and repeat until you reach the bottom.

Result: Row echelon form (REF) — a staircase pattern of zeros below the pivots.

Back substitution: With the triangular system in hand, solve for the last unknown first, then substitute back up.

Reduced row echelon form (RREF): Continue elimination upward to zero out entries above the pivots too, and scale each pivot to 1. The result directly reveals the solution (or the free variables).

The process in algebra: LU decomposition. Gaussian elimination factors $A = LU$ where $L$ is lower triangular (the elimination multipliers) and $U$ is upper triangular (the result after elimination). This is the fundamental factorization of linear algebra and is what numpy.linalg.solve actually computes under the hood.

Complexity: Gaussian elimination on an $n \times n$ system costs $O(n^3)$ operations. For sparse systems (most entries zero, common in physics and graphs), smarter algorithms reduce this dramatically.

Numerical stability: In finite-precision arithmetic, dividing by a tiny pivot amplifies errors explosively. Partial pivoting (always swap to bring the largest available entry to the pivot position) is the standard fix and is what all production libraries implement.

import numpy as np

def gaussian_elimination(A, b):
    """Solve Ax = b by Gaussian elimination with partial pivoting."""
    n = len(b)
    # Augmented matrix
    M = np.hstack([A.astype(float), b.reshape(-1, 1)])
    
    for col in range(n):
        # Partial pivoting: swap in the largest row below
        max_row = col + np.argmax(abs(M[col:, col]))
        M[[col, max_row]] = M[[max_row, col]]
        
        pivot = M[col, col]
        if abs(pivot) < 1e-12:
            raise ValueError("Matrix is singular")
        
        # Eliminate below
        for row in range(col + 1, n):
            factor = M[row, col] / pivot
            M[row, col:] -= factor * M[col, col:]
    
    # Back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (M[i, -1] - M[i, i+1:n] @ x[i+1:n]) / M[i, i]
    return x

A = np.array([[2., 1., -1.], [-3., -1., 2.], [-2., 1., 2.]])
b = np.array([8., -11., -3.])
print(gaussian_elimination(A, b))  # [2.  3. -1.]

Rank: The rank of matrix $A$ is the number of pivots in its row echelon form, which equals both the dimension of the column space and the dimension of the row space. Rank measures how much independent information is in the matrix.

Full rank: A square $n \times n$ matrix has full rank if $\text{rank}(A) = n$. A full-rank square matrix is invertible (non-singular): there exists a matrix $A^{-1}$ such that $AA^{-1} = A^{-1}A = I$.

Equivalent conditions for an $n \times n$ matrix $A$ to be invertible (all equivalent — if one holds, all hold):

• $\det(A) \neq 0$
• $A$ has $n$ pivots
• $A\mathbf{x} = \mathbf{0}$ has only the trivial solution $\mathbf{x} = \mathbf{0}$
• The columns of $A$ are linearly independent
• The rows of $A$ are linearly independent
• $A\mathbf{x} = \mathbf{b}$ has a unique solution for every b

The determinant $\det(A)$: For a $2 \times 2$ matrix, $\det\begin{pmatrix}a & b \ c & d\end{pmatrix} = ad - bc$. Geometrically, $|\det(A)|$ is the area scaling factor of the transformation represented by $A$: it tells you how much the matrix stretches (or squishes) areas. If $\det(A) = 0$, the transformation collapses at least one dimension — the matrix is singular. If $\det(A) < 0$, the transformation includes a reflection (orientation is flipped).

For $3 \times 3$, the determinant measures volume scaling. This geometric interpretation — determinant as the signed volume of the parallelepiped formed by the column vectors — is the cleanest way to understand why $\det(AB) = \det(A)\det(B)$: applying $B$ scales volumes by $\det(B)$, then $A$ scales by $\det(A)$, so total scaling is the product.

Practical note: Never invert a matrix to solve $A\mathbf{x} = \mathbf{b}$. Always use LU factorization (i.e., np.linalg.solve). Computing $A^{-1}$ explicitly is expensive and numerically inferior.

import numpy as np

A = np.array([[3., 1.], [2., 4.]])

print("det:", np.linalg.det(A))  # 3*4 - 1*2 = 10
print("rank:", np.linalg.matrix_rank(A))  # 2 (full rank)

A_inv = np.linalg.inv(A)
print("A @ A_inv:")   
print(np.round(A @ A_inv, 10))  # Identity matrix

# Singular matrix: columns are dependent
B = np.array([[1., 2.], [2., 4.]])  # row2 = 2*row1
print("det(B):", np.linalg.det(B))    # ~0
print("rank(B):", np.linalg.matrix_rank(B))  # 1