03. Linear Algebra

Previous 02. Complex Numbers
๐Ÿ“‹ Topic list
Next 04. Partial Differentiation

03. Linear Algebra

Boas Chapter 3 โ€” In physical sciences, linear algebra forms the foundation for nearly every field, including inertia tensors in mechanics, matrix formalism in quantum mechanics, and coupled oscillation analysis.

Learning Objectives

After completing this lesson, you will be able to:

  • Perform matrix operations (addition, multiplication, transpose, inverse) and calculate determinants
  • Solve systems of linear equations using Gaussian elimination and Cramer's rule
  • Find eigenvalues/eigenvectors and perform matrix diagonalization
  • Understand the spectral theorem for symmetric and Hermitian matrices, and utilize properties of orthogonal/unitary matrices
  • Determine the definiteness (positive/negative definite) of quadratic forms
  • Physics applications: Handle principal axis transformation of inertia tensors, normal modes of coupled oscillations, and matrix mechanics in quantum mechanics

1. Matrix Fundamentals

1.1 Matrices and Basic Operations

A matrix is a rectangular array of numbers. For an $m \times n$ matrix $A$, elements are denoted as $a_{ij}$:

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

Matrix addition: Element-wise addition for matrices of the same size

$$ (A + B)_{ij} = a_{ij} + b_{ij} $$

Scalar multiplication: Multiply all elements by a scalar

$$ (cA)_{ij} = c \cdot a_{ij} $$

Matrix multiplication: If $A$ is $m \times n$ and $B$ is $n \times p$, the result is $m \times p$

$$ (AB)_{ij} = \sum_{k=1}^{n} a_{ik} b_{kj} $$

Note: Matrix multiplication generally does not commute: $AB \neq BA$

import numpy as np

# ํ–‰๋ ฌ ๊ธฐ๋ณธ ์—ฐ์‚ฐ
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

print("A =\n", A)
print("B =\n", B)

# ํ–‰๋ ฌ ๊ณฑ์…ˆ
print("\nAB =\n", A @ B)
print("BA =\n", B @ A)
print("AB โ‰  BA:", not np.allclose(A @ B, B @ A))

# ํ–‰๋ ฌ ๊ณฑ์˜ ์„ฑ์งˆ
C = np.array([[1, 0], [2, 3]])
print("\n(AB)C =\n", (A @ B) @ C)
print("A(BC) =\n", A @ (B @ C))
print("๊ฒฐํ•ฉ๋ฒ•์น™ ์„ฑ๋ฆฝ:", np.allclose((A @ B) @ C, A @ (B @ C)))

1.2 Transpose and Conjugate Transpose

Transpose $A^T$: Exchange rows and columns

$$ (A^T)_{ij} = a_{ji} $$

Properties: - $(AB)^T = B^T A^T$ (order reversed) - $(A^T)^T = A$ - $(A + B)^T = A^T + B^T$

Conjugate transpose (adjoint) $A^\dagger$: Transpose + complex conjugate

$$ (A^\dagger)_{ij} = \overline{a_{ji}} $$

# ์ „์น˜ ํ–‰๋ ฌ
A = np.array([[1, 2, 3], [4, 5, 6]])
print("A =\n", A)
print("A^T =\n", A.T)

# ๋ณต์†Œ ํ–‰๋ ฌ์˜ ์ผค๋ ˆ ์ „์น˜
Z = np.array([[1+2j, 3-1j], [4j, 2+3j]])
print("\nZ =\n", Z)
print("Zโ€  =\n", Z.conj().T)

# (AB)^T = B^T A^T ํ™•์ธ
B = np.array([[1, 2], [3, 4], [5, 6]])
print("\n(AB)^T =\n", (A @ B).T)
print("B^T A^T =\n", B.T @ A.T)

1.3 Special Matrices

Identity matrix $I$: Diagonal elements are all 1

$$ I_{ij} = \delta_{ij} = \begin{cases} 1 & (i = j) \\ 0 & (i \neq j) \end{cases} $$

Diagonal matrix: Matrix with only diagonal elements non-zero

Symmetric matrix: $A^T = A$, i.e., $a_{ij} = a_{ji}$

Antisymmetric matrix: $A^T = -A$

Hermitian matrix: $A^\dagger = A$ (complex extension of symmetric matrices)

Orthogonal matrix: $A^T A = AA^T = I$, i.e., $A^{-1} = A^T$

Unitary matrix: $A^\dagger A = AA^\dagger = I$ (complex extension of orthogonal matrices)

# ์ง๊ต ํ–‰๋ ฌ ์˜ˆ: 2D ํšŒ์ „ ํ–‰๋ ฌ
theta = np.pi / 4  # 45๋„ ํšŒ์ „
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])

print("R =\n", R)
print("R^T R =\n", np.round(R.T @ R, 10))
print("det(R) =", np.linalg.det(R))  # det = +1 (proper rotation)

# ์—๋ฅด๋ฏธํŠธ ํ–‰๋ ฌ ์˜ˆ
H = np.array([[2, 1-1j], [1+1j, 3]])
print("\nH =\n", H)
print("Hโ€  =\n", H.conj().T)
print("์—๋ฅด๋ฏธํŠธ:", np.allclose(H, H.conj().T))

2. Determinants

2.1 Definition and Basic Properties

The determinant $\det(A)$ of an $n \times n$ square matrix $A$ is a scalar value.

2ร—2 determinant:

$$ \det\begin{pmatrix} a & b \\ c & d \end{pmatrix} = ad - bc $$

3ร—3 determinant (Sarrus's rule or cofactor expansion):

$$ \det\begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} = a_{11}(a_{22}a_{33} - a_{23}a_{32}) - a_{12}(a_{21}a_{33} - a_{23}a_{31}) + a_{13}(a_{21}a_{32} - a_{22}a_{31}) $$

Key properties:

  1. $\det(AB) = \det(A) \cdot \det(B)$
  2. $\det(A^T) = \det(A)$
  3. $\det(cA) = c^n \det(A)$ ($n \times n$ matrix)
  4. Exchanging rows (or columns) changes the sign
  5. If one row is a constant multiple of another, $\det = 0$
  6. $\det(A^{-1}) = 1/\det(A)$

2.2 Cofactor Expansion

Minor $M_{ij}$: Determinant of the $(n-1) \times (n-1)$ matrix obtained by deleting the $i$-th row and $j$-th column

Cofactor: $C_{ij} = (-1)^{i+j} M_{ij}$

Cofactor expansion of determinant (expanding along the $i$-th row):

$$ \det(A) = \sum_{j=1}^{n} a_{ij} C_{ij} = \sum_{j=1}^{n} (-1)^{i+j} a_{ij} M_{ij} $$

import numpy as np
from numpy.linalg import det

# ํ–‰๋ ฌ์‹ ๊ณ„์‚ฐ
A = np.array([[2, 1, 3],
              [0, -1, 2],
              [4, 3, 1]])

print("det(A) =", det(A))

# ์—ฌ์ธ์ˆ˜ ์ „๊ฐœ๋ฅผ ์ˆ˜๋™์œผ๋กœ ๊ณ„์‚ฐ (1ํ–‰ ๊ธฐ์ค€)
# det = 2*(-1-6) - 1*(0-8) + 3*(0+4)
manual = 2*(-1*1 - 2*3) - 1*(0*1 - 2*4) + 3*(0*3 - (-1)*4)
print("์ˆ˜๋™ ๊ณ„์‚ฐ:", manual)

# ํ–‰๋ ฌ์‹์˜ ์„ฑ์งˆ ํ™•์ธ
B = np.array([[1, 2, 0],
              [3, 1, -1],
              [2, 0, 4]])

print(f"\ndet(A) = {det(A):.4f}")
print(f"det(B) = {det(B):.4f}")
print(f"det(AB) = {det(A @ B):.4f}")
print(f"det(A)*det(B) = {det(A)*det(B):.4f}")

2.3 Geometric Interpretation of Determinants

The absolute value $|\det(A)|$ of the determinant of an $n \times n$ matrix $A$ is the volume of the parallelepiped formed by the row vectors (or column vectors).

  • 2D: $|\det(A)|$ = area of the parallelogram formed by two vectors
  • 3D: $|\det(A)|$ = volume of the parallelepiped formed by three vectors

The sign of $\det(A)$ indicates orientation: - $\det > 0$: Right-handed coordinate system preserved - $\det < 0$: Coordinate system inverted (includes mirror reflection)

import matplotlib.pyplot as plt

# 2D์—์„œ์˜ ํ–‰๋ ฌ์‹ = ๋„“์ด
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ๋ณ€ํ™˜ ์ „: ๋‹จ์œ„ ์ •์‚ฌ๊ฐํ˜•
square = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]).T

A = np.array([[2, 1], [0, 3]])  # det = 6

# ๋ณ€ํ™˜ ํ›„
transformed = A @ square

for ax, shape, title in [(axes[0], square, '๋‹จ์œ„ ์ •์‚ฌ๊ฐํ˜• (๋„“์ด=1)'),
                          (axes[1], transformed, f'๋ณ€ํ™˜ ํ›„ (๋„“์ด=|det|={abs(det(A)):.0f})')]:
    ax.fill(shape[0], shape[1], alpha=0.3, color='blue')
    ax.plot(shape[0], shape[1], 'b-', linewidth=2)
    ax.set_xlim(-1, 8)
    ax.set_ylim(-1, 5)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    ax.set_title(title)
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.axvline(x=0, color='k', linewidth=0.5)

plt.tight_layout()
plt.savefig('determinant_area.png', dpi=100, bbox_inches='tight')
plt.close()
print(f"det(A) = {det(A):.0f}: ๋„“์ด๊ฐ€ {abs(det(A)):.0f}๋ฐฐ ํ™•๋Œ€")

3. Inverse Matrices and Systems of Linear Equations

3.1 Inverse Matrix

The inverse matrix $A^{-1}$ of a square matrix $A$ satisfies:

$$ AA^{-1} = A^{-1}A = I $$

Necessary and sufficient condition for the inverse to exist: $\det(A) \neq 0$ (non-singular matrix)

Inverse using cofactors:

$$ A^{-1} = \frac{1}{\det(A)} \text{adj}(A) $$

where $\text{adj}(A)$ is the adjugate matrix, the transpose of the cofactor matrix: $\text{adj}(A)_{ij} = C_{ji}$

2ร—2 inverse:

$$ \begin{pmatrix} a & b \\ c & d \end{pmatrix}^{-1} = \frac{1}{ad - bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} $$

A = np.array([[2, 1], [5, 3]])
A_inv = np.linalg.inv(A)

print("A =\n", A)
print("A^{-1} =\n", A_inv)
print("A @ A^{-1} =\n", np.round(A @ A_inv, 10))

# ์ˆ˜๋ฐ˜ ํ–‰๋ ฌ์„ ์ด์šฉํ•œ ์ˆ˜๋™ ๊ณ„์‚ฐ
d = det(A)  # 2*3 - 1*5 = 1
adj_A = np.array([[3, -1], [-5, 2]])  # ์—ฌ์ธ์ˆ˜ ์ „์น˜
manual_inv = adj_A / d
print("\n์ˆ˜๋™ ๊ณ„์‚ฐ:\n", manual_inv)

3.2 Gaussian Elimination

Solve the system of linear equations $A\mathbf{x} = \mathbf{b}$ by performing row operations on the augmented matrix $(A | \mathbf{b})$.

Allowed elementary row operations: 1. Exchange two rows 2. Multiply a row by a non-zero constant 3. Add a constant multiple of one row to another row

Goal: Transform to upper triangular form (or reduced row echelon form) followed by back substitution

def gauss_eliminate(A_aug, verbose=True):
    """๊ฐ€์šฐ์Šค ์†Œ๊ฑฐ๋ฒ•์œผ๋กœ ์—ฐ๋ฆฝ๋ฐฉ์ •์‹์„ ํ’‰๋‹ˆ๋‹ค.

    A_aug: ํ™•๋Œ€ ํ–‰๋ ฌ [A | b]
    """
    A = A_aug.astype(float).copy()
    n = A.shape[0]

    # ์ „์ง„ ์†Œ๊ฑฐ (forward elimination)
    for col in range(n):
        # ํ”ผ๋ด‡ ์„ ํƒ (๋ถ€๋ถ„ ํ”ผ๋ด‡ํŒ…)
        max_row = col + np.argmax(np.abs(A[col:, col]))
        if max_row != col:
            A[[col, max_row]] = A[[max_row, col]]

        if abs(A[col, col]) < 1e-12:
            print(f"๊ฒฝ๊ณ : ํ”ผ๋ด‡์ด 0์— ๊ฐ€๊น์Šต๋‹ˆ๋‹ค (์—ด {col})")
            continue

        # ์•„๋ž˜ ํ–‰๋“ค ์†Œ๊ฑฐ
        for row in range(col + 1, n):
            factor = A[row, col] / A[col, col]
            A[row] -= factor * A[col]

    if verbose:
        print("์ƒ์‚ผ๊ฐ ํ–‰๋ ฌ:\n", np.round(A, 4))

    # ํ›„์ง„ ๋Œ€์ž… (back substitution)
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (A[i, -1] - A[i, i+1:n] @ x[i+1:n]) / A[i, i]

    return x

# ์˜ˆ์ œ: 3x3 ์—ฐ๋ฆฝ๋ฐฉ์ •์‹
# 2x + y - z = 8
# -3x - y + 2z = -11
# -2x + y + 2z = -3
A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])
b = np.array([8, -11, -3])

A_aug = np.column_stack([A, b])
x = gauss_eliminate(A_aug)
print(f"\nํ•ด: x = {x}")
print(f"๊ฒ€์ฆ Ax = {A @ x}")

# NumPy ๋‚ด์žฅ ํ•จ์ˆ˜์™€ ๋น„๊ต
x_np = np.linalg.solve(A, b)
print(f"NumPy solve: {x_np}")

3.3 Cramer's Rule

For a system of linear equations $A\mathbf{x} = \mathbf{b}$, if $\det(A) \neq 0$, each unknown is:

$$ x_i = \frac{\det(A_i)}{\det(A)} $$

where $A_i$ is the matrix obtained by replacing the $i$-th column of $A$ with $\mathbf{b}$.

Note: Cramer's rule is theoretically important but computationally less efficient than Gaussian elimination. For large $n$, Gaussian elimination ($O(n^3)$) is much faster than Cramer's rule ($O(n \cdot n!)$).

def cramer(A, b):
    """ํฌ๋ž˜๋จธ ๋ฒ•์น™์œผ๋กœ Ax = b๋ฅผ ํ’‰๋‹ˆ๋‹ค."""
    n = len(b)
    d = det(A)
    if abs(d) < 1e-12:
        raise ValueError("ํ–‰๋ ฌ์‹์ด 0: ์œ ์ผํ•œ ํ•ด๊ฐ€ ์กด์žฌํ•˜์ง€ ์•Š์Œ")

    x = np.zeros(n)
    for i in range(n):
        A_i = A.copy()
        A_i[:, i] = b
        x[i] = det(A_i) / d
    return x

# ์˜ˆ์ œ
A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])
b = np.array([8, -11, -3])

x = cramer(A, b)
print(f"ํฌ๋ž˜๋จธ ๋ฒ•์น™ ํ•ด: {x}")

3.4 Rank of a Matrix

The rank of a matrix is the maximum number of linearly independent rows (or columns).

$$ \text{rank}(A) = \text{dim}(\text{Col}(A)) = \text{dim}(\text{Row}(A)) $$

Existence of solutions for systems of linear equations: - $\text{rank}(A) = \text{rank}(A|b) = n$: Unique solution - $\text{rank}(A) = \text{rank}(A|b) < n$: Infinitely many solutions - $\text{rank}(A) < \text{rank}(A|b)$: No solution (inconsistent)

# ๊ณ„์ˆ˜ ์˜ˆ์ œ
A_full = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])  # rank 3
A_deficient = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])  # rank 2 (3ํ–‰ = ํ–‰1 + 2*ํ–‰2 ์•„๋‹˜, ํ•˜์ง€๋งŒ 1+3=2*2)

print(f"rank(A_full) = {np.linalg.matrix_rank(A_full)}")
print(f"rank(A_deficient) = {np.linalg.matrix_rank(A_deficient)}")
print(f"det(A_full) = {det(A_full):.4f}")
print(f"det(A_deficient) = {det(A_deficient):.4f}")

4. Eigenvalues and Eigenvectors

4.1 Definition and Characteristic Equation

For a square matrix $A$, find a non-zero vector $\mathbf{v}$ and scalar $\lambda$ that satisfy:

$$ A\mathbf{v} = \lambda\mathbf{v} $$

  • $\lambda$: eigenvalue
  • $\mathbf{v}$: eigenvector

Geometric interpretation: Under the linear transformation by $A$, an eigenvector is a special vector whose direction remains unchanged. The eigenvalue is the scaling factor in that direction.

Characteristic equation:

$$ \det(A - \lambda I) = 0 $$

This is an $n$-th degree polynomial in $\lambda$, called the characteristic polynomial.

import numpy as np
from numpy.linalg import eig

# 2x2 ํ–‰๋ ฌ์˜ ๊ณ ์œ ๊ฐ’/๊ณ ์œ ๋ฒกํ„ฐ
A = np.array([[4, 2], [1, 3]])

eigenvalues, eigenvectors = eig(A)
print("๊ณ ์œ ๊ฐ’:", eigenvalues)
print("๊ณ ์œ ๋ฒกํ„ฐ (์—ด๋ฒกํ„ฐ):\n", eigenvectors)

# ๊ฒ€์ฆ: Av = ฮปv
for i in range(len(eigenvalues)):
    lam = eigenvalues[i]
    v = eigenvectors[:, i]
    Av = A @ v
    lam_v = lam * v
    print(f"\nฮป_{i+1} = {lam:.4f}")
    print(f"v_{i+1} = {v}")
    print(f"Av = {Av}")
    print(f"ฮปv = {lam_v}")
    print(f"Av โ‰ˆ ฮปv: {np.allclose(Av, lam_v)}")

4.2 Manual Calculation Example

Find eigenvalues/eigenvectors of $A = \begin{pmatrix} 3 & 1 \\ 0 & 2 \end{pmatrix}$:

Characteristic equation: $$ \det(A - \lambda I) = \det\begin{pmatrix} 3-\lambda & 1 \\ 0 & 2-\lambda \end{pmatrix} = (3-\lambda)(2-\lambda) = 0 $$

$\lambda_1 = 3$, $\lambda_2 = 2$

Eigenvector ($\lambda_1 = 3$): $$ (A - 3I)\mathbf{v} = \begin{pmatrix} 0 & 1 \\ 0 & -1 \end{pmatrix}\mathbf{v} = 0 \implies v_2 = 0 \implies \mathbf{v}_1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix} $$

Eigenvector ($\lambda_2 = 2$): $$ (A - 2I)\mathbf{v} = \begin{pmatrix} 1 & 1 \\ 0 & 0 \end{pmatrix}\mathbf{v} = 0 \implies v_1 + v_2 = 0 \implies \mathbf{v}_2 = \begin{pmatrix} -1 \\ 1 \end{pmatrix} $$

import sympy as sp

# SymPy๋ฅผ ์ด์šฉํ•œ ํ•ด์„์  ๊ณ„์‚ฐ
A_sym = sp.Matrix([[3, 1], [0, 2]])
print("ํŠน์„ฑ ๋‹คํ•ญ์‹:", A_sym.charpoly().as_expr())
print("๊ณ ์œ ๊ฐ’:", A_sym.eigenvals())  # {3: 1, 2: 1} (๊ณ ์œ ๊ฐ’: ์ค‘๋ณต๋„)
print("๊ณ ์œ ๋ฒกํ„ฐ:", A_sym.eigenvects())

4.3 Matrix Diagonalization

An $n \times n$ matrix $A$ is diagonalizable if it has $n$ linearly independent eigenvectors.

Construct a matrix $P = [\mathbf{v}_1 | \mathbf{v}_2 | \cdots | \mathbf{v}_n]$ with eigenvectors as columns:

$$ P^{-1}AP = D = \text{diag}(\lambda_1, \lambda_2, \ldots, \lambda_n) $$

Usefulness of diagonalization:

  • Powers: $A^k = P D^k P^{-1}$, where $D^k$ is computed by raising diagonal elements to the $k$-th power
  • Exponential: $e^{At} = P \, e^{Dt} \, P^{-1}$, where $e^{Dt} = \text{diag}(e^{\lambda_1 t}, \ldots, e^{\lambda_n t})$
  • Coupled ODEs: Decouple $\dot{\mathbf{x}} = A\mathbf{x}$ into independent scalar ODEs
# ๋Œ€๊ฐํ™” ์˜ˆ์ œ
A = np.array([[4, 2], [1, 3]])
eigenvalues, P = eig(A)
D = np.diag(eigenvalues)

print("P (๊ณ ์œ ๋ฒกํ„ฐ ํ–‰๋ ฌ):\n", P)
print("D (๋Œ€๊ฐ ํ–‰๋ ฌ):\n", D)

# ๊ฒ€์ฆ: A = P D P^{-1}
P_inv = np.linalg.inv(P)
A_reconstructed = P @ D @ P_inv
print("\nP D P^{-1} =\n", np.round(A_reconstructed, 10))
print("A์™€ ์ผ์น˜:", np.allclose(A, A_reconstructed))

# ์‘์šฉ: A^10 ๊ณ„์‚ฐ
A_power_10 = P @ np.diag(eigenvalues**10) @ P_inv
print(f"\nA^10 =\n{np.round(A_power_10.real, 2)}")
print(f"์ง์ ‘ ๊ณ„์‚ฐ:\n{np.round(np.linalg.matrix_power(A, 10).astype(float), 2)}")

4.4 Properties of Symmetric and Hermitian Matrices

Spectral Theorem:

  1. Real symmetric matrices ($A^T = A$):
  2. All eigenvalues are real
  3. Eigenvectors corresponding to distinct eigenvalues are orthogonal

  4. Diagonalizable by an orthogonal matrix $Q$: $Q^T A Q = D$, $Q^T = Q^{-1}$

  5. Hermitian matrices ($A^\dagger = A$):

  6. All eigenvalues are real
  7. Eigenvectors corresponding to distinct eigenvalues are orthogonal
  8. Diagonalizable by a unitary matrix $U$: $U^\dagger A U = D$, $U^\dagger = U^{-1}$
# ์‹ค๋Œ€์นญ ํ–‰๋ ฌ์˜ ์„ฑ์งˆ
A_sym = np.array([[4, 1, 2],
                   [1, 3, 1],
                   [2, 1, 5]])

eigenvalues, eigenvectors = eig(A_sym)
print("๋Œ€์นญ ํ–‰๋ ฌ์˜ ๊ณ ์œ ๊ฐ’:", eigenvalues)
print("๋ชจ๋‘ ์‹ค์ˆ˜:", np.allclose(eigenvalues.imag, 0))

# ๊ณ ์œ ๋ฒกํ„ฐ์˜ ์ง๊ต์„ฑ ํ™•์ธ
print("\n๊ณ ์œ ๋ฒกํ„ฐ ๋‚ด์ :")
for i in range(3):
    for j in range(i+1, 3):
        dot = eigenvectors[:, i] @ eigenvectors[:, j]
        print(f"  v{i+1} ยท v{j+1} = {dot:.6e}")

# ์ง๊ต ํ–‰๋ ฌ๋กœ ๋Œ€๊ฐํ™”
Q = eigenvectors
print(f"\nQ^T Q =\n{np.round(Q.T @ Q, 10)}")
print("Q๋Š” ์ง๊ต ํ–‰๋ ฌ:", np.allclose(Q.T @ Q, np.eye(3)))

4.5 Summary of Major Matrix Decompositions

Decomposition Form Condition Physics Application
Eigenvalue decomposition $A = PDP^{-1}$ $n$ independent eigenvectors Oscillation modes, stability
Orthogonal diagonalization $A = QDQ^T$ Symmetric matrix Principal axis theorem, inertia tensor
SVD $A = U\Sigma V^T$ Any $m \times n$ Data reduction, least squares
LU decomposition $A = LU$ Square matrix Efficient solution of linear systems
Cholesky $A = LL^T$ Positive definite symmetric Statistics, optimization

5. Quadratic Forms and Definiteness

5.1 Definition of Quadratic Forms

A quadratic form for a vector $\mathbf{x} \in \mathbb{R}^n$:

$$ Q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x} = \sum_{i,j} a_{ij} x_i x_j $$

where $A$ can always be treated as a symmetric matrix (the antisymmetric part does not contribute to the quadratic form, so replace $A \to (A + A^T)/2$).

2-variable example:

$$ Q(x, y) = \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} a & b \\ b & c \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = ax^2 + 2bxy + cy^2 $$

5.2 Determining Definiteness

Definiteness of a symmetric matrix $A$:

Classification Condition Eigenvalue Condition
Positive definite $\mathbf{x}^T A \mathbf{x} > 0$ (for all $\mathbf{x} \neq 0$) All $\lambda_i > 0$
Positive semidefinite $\mathbf{x}^T A \mathbf{x} \geq 0$ All $\lambda_i \geq 0$
Negative definite $\mathbf{x}^T A \mathbf{x} < 0$ (for all $\mathbf{x} \neq 0$) All $\lambda_i < 0$
Indefinite Can be both positive and negative Mixed positive/negative eigenvalues

Sylvester's criterion: A symmetric matrix is positive definite if and only if all leading principal minors are positive:

$$ a_{11} > 0, \quad \det\begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} > 0, \quad \det\begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} > 0, \quad \ldots $$

def check_definiteness(A):
    """๋Œ€์นญ ํ–‰๋ ฌ์˜ ์ •์น˜์„ฑ์„ ํŒ๋ณ„ํ•ฉ๋‹ˆ๋‹ค."""
    eigenvalues = np.linalg.eigvalsh(A)  # ๋Œ€์นญ ํ–‰๋ ฌ ์ „์šฉ
    print(f"๊ณ ์œ ๊ฐ’: {eigenvalues}")

    if all(eigenvalues > 0):
        return "์–‘์ •์น˜ (positive definite)"
    elif all(eigenvalues >= 0):
        return "์–‘๋ฐ˜์ •์น˜ (positive semidefinite)"
    elif all(eigenvalues < 0):
        return "์Œ์ •์น˜ (negative definite)"
    elif all(eigenvalues <= 0):
        return "์Œ๋ฐ˜์ •์น˜ (negative semidefinite)"
    else:
        return "๋ถ€์ •์น˜ (indefinite)"

# ์˜ˆ์ œ
matrices = {
    "์–‘์ •์น˜": np.array([[2, -1], [-1, 2]]),
    "๋ถ€์ •์น˜": np.array([[1, 3], [3, 1]]),
    "์–‘๋ฐ˜์ •์น˜": np.array([[1, 1], [1, 1]]),
}

for name, M in matrices.items():
    result = check_definiteness(M)
    print(f"{name} ํ–‰๋ ฌ: {result}\n")

5.3 Principal Axis Theorem

When a quadratic form is transformed to principal axis coordinates, cross terms disappear.

If a symmetric matrix $A = Q D Q^T$:

$$ Q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x} = \mathbf{y}^T D \mathbf{y} = \lambda_1 y_1^2 + \lambda_2 y_2^2 + \cdots + \lambda_n y_n^2 $$

where $\mathbf{y} = Q^T \mathbf{x}$ (principal axis coordinates).

This is directly applied in physics to principal axis transformation of inertia tensors, principal stress directions of stress tensors, etc.

import matplotlib.pyplot as plt

# ์ด์ฐจ๊ณก๋ฉด๊ณผ ์ฃผ์ถ• ๋ณ€ํ™˜ ์‹œ๊ฐํ™”
A = np.array([[5, 2], [2, 2]])
eigenvalues, Q = np.linalg.eigh(A)

print(f"์›๋ž˜ ์ด์ฐจํ˜•์‹: 5xยฒ + 4xy + 2yยฒ")
print(f"๊ณ ์œ ๊ฐ’: ฮปโ‚={eigenvalues[0]:.2f}, ฮปโ‚‚={eigenvalues[1]:.2f}")
print(f"์ฃผ์ถ• ์ขŒํ‘œ: {eigenvalues[0]:.2f}yโ‚ยฒ + {eigenvalues[1]:.2f}yโ‚‚ยฒ")

# ํƒ€์› ์‹œ๊ฐํ™”
theta = np.linspace(0, 2*np.pi, 200)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ์›๋ž˜ ์ขŒํ‘œ๊ณ„์—์„œ์˜ ํƒ€์› (5xยฒ + 4xy + 2yยฒ = 1)
# ๋งค๊ฐœ๋ณ€์ˆ˜ ํ‘œํ˜„ ์ด์šฉ
t = np.linspace(0, 2*np.pi, 200)
# ์ฃผ์ถ• ์ขŒํ‘œ์—์„œ ํƒ€์›: ฮปโ‚yโ‚ยฒ + ฮปโ‚‚yโ‚‚ยฒ = 1
y1 = np.cos(t) / np.sqrt(eigenvalues[0])
y2 = np.sin(t) / np.sqrt(eigenvalues[1])
# ์›๋ž˜ ์ขŒํ‘œ๋กœ ๋ณ€ํ™˜: x = Q y
xy = Q @ np.array([y1, y2])

axes[0].plot(xy[0], xy[1], 'b-', linewidth=2)
# ์ฃผ์ถ• ๋ฐฉํ–ฅ
for i in range(2):
    v = Q[:, i] / np.sqrt(eigenvalues[i])
    axes[0].arrow(0, 0, v[0], v[1], head_width=0.03, color=['red', 'green'][i],
                  linewidth=2, label=f'์ฃผ์ถ• {i+1} (ฮป={eigenvalues[i]:.2f})')
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)
axes[0].legend()
axes[0].set_title('์›๋ž˜ ์ขŒํ‘œ๊ณ„: 5xยฒ + 4xy + 2yยฒ = 1')

# ์ฃผ์ถ• ์ขŒํ‘œ์—์„œ์˜ ํƒ€์› (์ถ• ์ •๋ ฌ)
y1 = np.cos(t) / np.sqrt(eigenvalues[0])
y2 = np.sin(t) / np.sqrt(eigenvalues[1])
axes[1].plot(y1, y2, 'b-', linewidth=2)
axes[1].axhline(y=0, color='green', linewidth=1.5, label=f'yโ‚‚์ถ• (ฮป={eigenvalues[1]:.2f})')
axes[1].axvline(x=0, color='red', linewidth=1.5, label=f'yโ‚์ถ• (ฮป={eigenvalues[0]:.2f})')
axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)
axes[1].legend()
axes[1].set_title(f'์ฃผ์ถ• ์ขŒํ‘œ: {eigenvalues[0]:.2f}yโ‚ยฒ + {eigenvalues[1]:.2f}yโ‚‚ยฒ = 1')

plt.tight_layout()
plt.savefig('principal_axes.png', dpi=100, bbox_inches='tight')
plt.close()

6. Physics Applications

6.1 Moment of Inertia Tensor and Principal Axes

The moment of inertia tensor of a rigid body is a $3 \times 3$ symmetric matrix:

$$ I = \begin{pmatrix} I_{xx} & -I_{xy} & -I_{xz} \\ -I_{yx} & I_{yy} & -I_{yz} \\ -I_{zx} & -I_{zy} & I_{zz} \end{pmatrix} $$

where: - $I_{xx} = \sum m_i (y_i^2 + z_i^2)$ (moment of inertia) - $I_{xy} = \sum m_i x_i y_i$ (product of inertia)

Principal axis transformation: Diagonalizing $I$ yields the principal moments of inertia $I_1, I_2, I_3$ and principal axes.

$$ Q^T I Q = \text{diag}(I_1, I_2, I_3) $$

In the principal axis coordinate system, angular momentum simplifies to $L_i = I_i \omega_i$.

# ๊ด€์„ฑ ํ…์„œ ์˜ˆ์ œ: ์ •์œก๋ฉด์ฒด ๊ผญ์ง“์ ์— ์งˆ๋Ÿ‰์ด ์žˆ๋Š” ๊ฒฝ์šฐ
# ์ •์œก๋ฉด์ฒด ๊ผญ์ง“์  ์ขŒํ‘œ (๋ณ€์˜ ๊ธธ์ด 2, ์ค‘์‹ฌ ์›์ )
vertices = np.array([
    [ 1,  1,  1], [ 1,  1, -1], [ 1, -1,  1], [ 1, -1, -1],
    [-1,  1,  1], [-1,  1, -1], [-1, -1,  1], [-1, -1, -1]
], dtype=float)

m = 1.0  # ๊ฐ ๊ผญ์ง“์ ์˜ ์งˆ๋Ÿ‰

# ๊ด€์„ฑ ํ…์„œ ๊ณ„์‚ฐ
I_tensor = np.zeros((3, 3))
for r in vertices:
    I_tensor += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r))

print("๊ด€์„ฑ ํ…์„œ:\n", I_tensor)

# ๋Œ€๊ฐํ™” โ†’ ์ฃผ๊ด€์„ฑ ๋ชจ๋ฉ˜ํŠธ
eigenvalues, eigenvectors = np.linalg.eigh(I_tensor)
print(f"\n์ฃผ๊ด€์„ฑ ๋ชจ๋ฉ˜ํŠธ: Iโ‚={eigenvalues[0]:.2f}, Iโ‚‚={eigenvalues[1]:.2f}, Iโ‚ƒ={eigenvalues[2]:.2f}")
print("์ฃผ์ถ• ๋ฐฉํ–ฅ:\n", eigenvectors)

6.2 Coupled Oscillations

A system with two masses connected by springs:

$$ M\ddot{\mathbf{x}} = -K\mathbf{x} $$

where $M$ is the mass matrix and $K$ is the stiffness matrix.

Normal frequencies and normal modes form a generalized eigenvalue problem:

$$ K\mathbf{v} = \omega^2 M\mathbf{v} $$

# ์—ฐ์„ฑ ์ง„๋™ ์˜ˆ์ œ: ๋‘ ์งˆ๋Ÿ‰-์Šคํ”„๋ง ๊ณ„
# mโ‚ = mโ‚‚ = m, ์Šคํ”„๋ง ์ƒ์ˆ˜ kโ‚ = kโ‚‚ = kโ‚ƒ = k
# ์šด๋™๋ฐฉ์ •์‹: m*xโ‚'' = -k*xโ‚ + k*(xโ‚‚-xโ‚) = -2k*xโ‚ + k*xโ‚‚
#             m*xโ‚‚'' = -k*(xโ‚‚-xโ‚) - k*xโ‚‚ = k*xโ‚ - 2k*xโ‚‚

m, k = 1.0, 1.0

M_mat = m * np.eye(2)
K_mat = np.array([[2*k, -k],
                   [-k, 2*k]])

# ์ผ๋ฐ˜ํ™” ๊ณ ์œ ๊ฐ’ ๋ฌธ์ œ: K v = ฯ‰ยฒ M v
# M = mI์ด๋ฏ€๋กœ (1/m)K v = ฯ‰ยฒ v
eigenvalues, eigenvectors = np.linalg.eigh(K_mat / m)
omega = np.sqrt(eigenvalues)

print("๊ฐ•์„ฑ ํ–‰๋ ฌ K:\n", K_mat)
print(f"\n๊ณ ์œ ์ง„๋™์ˆ˜: ฯ‰โ‚ = {omega[0]:.4f}, ฯ‰โ‚‚ = {omega[1]:.4f}")
print(f"์ฃผ๊ธฐ: Tโ‚ = {2*np.pi/omega[0]:.4f}, Tโ‚‚ = {2*np.pi/omega[1]:.4f}")
print(f"\n๊ณ ์œ ๋ชจ๋“œ 1 (๋™์ƒ): {eigenvectors[:, 0]}")
print(f"๊ณ ์œ ๋ชจ๋“œ 2 (์—ญ์ƒ): {eigenvectors[:, 1]}")

# ์‹œ๊ฐ„ ๋ฐœ์ „ ์‹œ๊ฐํ™”
t = np.linspace(0, 20, 500)
# ์ดˆ๊ธฐ ์กฐ๊ฑด: xโ‚(0) = 1, xโ‚‚(0) = 0, ์†๋„ 0
# ๋ชจ๋“œ ์ขŒํ‘œ๋กœ ๋ถ„ํ•ด
eta = eigenvectors.T @ np.array([1, 0])
x = np.zeros((2, len(t)))
for i in range(2):
    x += np.outer(eigenvectors[:, i], eta[i] * np.cos(omega[i] * t))

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(t, x[0], label='xโ‚(t)', linewidth=1.5)
ax.plot(t, x[1], label='xโ‚‚(t)', linewidth=1.5)
ax.set_xlabel('์‹œ๊ฐ„ t')
ax.set_ylabel('๋ณ€์œ„')
ax.set_title('์—ฐ์„ฑ ์ง„๋™: ๋งฅ๋†€์ด(beat) ํ˜„์ƒ')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('coupled_oscillations.png', dpi=100, bbox_inches='tight')
plt.close()

6.3 Matrices in Quantum Mechanics

In quantum mechanics, observables are represented as Hermitian operators, which in finite dimensions are Hermitian matrices.

Pauli matrices for spin-1/2 particles:

$$ \sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} $$

These are Hermitian ($\sigma_i^\dagger = \sigma_i$), and each has eigenvalues $\pm 1$.

# ํŒŒ์šธ๋ฆฌ ํ–‰๋ ฌ
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)

for name, sigma in [('ฯƒ_x', sigma_x), ('ฯƒ_y', sigma_y), ('ฯƒ_z', sigma_z)]:
    vals, vecs = np.linalg.eigh(sigma)
    print(f"{name}:")
    print(f"  ์—๋ฅด๋ฏธํŠธ: {np.allclose(sigma, sigma.conj().T)}")
    print(f"  ๊ณ ์œ ๊ฐ’: {vals}")
    print(f"  tr(ฯƒ) = {np.trace(sigma):.0f}, det(ฯƒ) = {np.linalg.det(sigma):.0f}")
    print()

# ๊ตํ™˜ ๊ด€๊ณ„: [ฯƒ_i, ฯƒ_j] = 2i ฮต_ijk ฯƒ_k
comm_xy = sigma_x @ sigma_y - sigma_y @ sigma_x
print("[ฯƒ_x, ฯƒ_y] = 2i ฯƒ_z:")
print(np.round(comm_xy, 10))
print("2i ฯƒ_z:")
print(2j * sigma_z)
print("์ผ์น˜:", np.allclose(comm_xy, 2j * sigma_z))

7. Linear Algebra with Python

7.1 Using NumPy and SciPy

import numpy as np
from scipy import linalg

# ===== ๊ธฐ๋ณธ ์—ฐ์‚ฐ =====
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])

# ํ–‰๋ ฌ์‹, ์—ญํ–‰๋ ฌ
print(f"det(A) = {np.linalg.det(A):.4f}")
print(f"A^(-1) =\n{np.linalg.inv(A)}")

# ๊ณ ์œ ๊ฐ’ ๋ถ„ํ•ด
vals, vecs = np.linalg.eig(A)
print(f"\n๊ณ ์œ ๊ฐ’: {vals}")

# ===== ๋‹ค์–‘ํ•œ ํ–‰๋ ฌ ๋ถ„ํ•ด =====

# LU ๋ถ„ํ•ด
P, L, U = linalg.lu(A)
print(f"\nLU ๋ถ„ํ•ด:")
print(f"P =\n{P}")
print(f"L =\n{np.round(L, 4)}")
print(f"U =\n{np.round(U, 4)}")
print(f"PLU = A: {np.allclose(P @ L @ U, A)}")

# SVD (Singular Value Decomposition)
U_svd, s, Vt = np.linalg.svd(A)
print(f"\nํŠน์ด๊ฐ’: {s}")

# QR ๋ถ„ํ•ด
Q, R = np.linalg.qr(A)
print(f"\nQR ๋ถ„ํ•ด: Q ์ง๊ต์„ฑ {np.allclose(Q.T @ Q, np.eye(3))}")

# ์—ฐ๋ฆฝ๋ฐฉ์ •์‹ ํ’€๊ธฐ
b = np.array([1, 2, 3])
x = np.linalg.solve(A, b)
print(f"\nAx = b์˜ ํ•ด: {x}")
print(f"๊ฒ€์ฆ: Ax = {A @ x}")

7.2 Symbolic Computation with SymPy

import sympy as sp

# ๊ธฐํ˜ธ ํ–‰๋ ฌ
a, b, c, d = sp.symbols('a b c d')
M = sp.Matrix([[a, b], [c, d]])

print("์ผ๋ฐ˜ 2x2 ํ–‰๋ ฌ:")
print(f"  det = {M.det()}")
print(f"  trace = {M.trace()}")
print(f"  ํŠน์„ฑ ๋‹คํ•ญ์‹: {M.charpoly().as_expr()}")

# ์ผ€์ผ๋ฆฌ-ํ•ด๋ฐ€ํ„ด ์ •๋ฆฌ: ๋ชจ๋“  ํ–‰๋ ฌ์€ ์ž์‹ ์˜ ํŠน์„ฑ ๋‹คํ•ญ์‹์„ ๋งŒ์กฑ
# Aยฒ - (a+d)A + (ad-bc)I = 0
lam = sp.Symbol('lambda')
char_poly = M.charpoly(lam)
print(f"\nํŠน์„ฑ ๋‹คํ•ญ์‹ p(ฮป) = {char_poly.as_expr()}")
print("์ผ€์ผ๋ฆฌ-ํ•ด๋ฐ€ํ„ด: p(M) = 0 ํ™•์ธ")
result = M**2 - (a+d)*M + (a*d - b*c)*sp.eye(2)
print(f"  p(M) = {sp.simplify(result)}")

Exercises

Basic Problems

1. Find the determinant and inverse of the following matrix: $$ A = \begin{pmatrix} 1 & 2 & 1 \\ 3 & 1 & -1 \\ 2 & 0 & 3 \end{pmatrix} $$

2. Solve the following system of equations using Gaussian elimination: $$ x + 2y + z = 4, \quad 3x + y - z = 2, \quad 2x + 3z = 7 $$

3. Find the eigenvalues and eigenvectors of the following matrix: $$ B = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} $$

Applied Problems

4. Determine whether the quadratic form $Q(x,y) = 3x^2 + 2xy + 3y^2$ is positive definite/negative definite/indefinite, and perform a principal axis transformation.

5. Two objects with masses $m_1 = 1$, $m_2 = 2$ are connected by three springs with spring constants $k_1 = 3$, $k_2 = 4$, $k_3 = 2$ (wallโ€“kโ‚โ€“mโ‚โ€“kโ‚‚โ€“mโ‚‚โ€“kโ‚ƒโ€“wall). Find the normal frequencies and normal modes.

6. For a spin-1/2 measurement operator in the magnetic field direction $\hat{n} = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)$, $\hat{S}_n = \frac{\hbar}{2}(\sigma_x \sin\theta\cos\phi + \sigma_y \sin\theta\sin\phi + \sigma_z \cos\theta)$, find the eigenvalues and eigenvectors.

References

  • Boas, Mathematical Methods in the Physical Sciences, Chapter 3
  • Strang, Introduction to Linear Algebra
  • Arfken & Weber, Mathematical Methods for Physicists, Chapters 3-4
  • Code for this lesson: See examples/Math_for_AI/01_vector_matrix_ops.py, 02_svd_pca.py

Next Lesson

04. Partial Differentiation covers calculus of multivariable functions.

Previous 02. Complex Numbers
๐Ÿ“‹ Topic list
Next 04. Partial Differentiation
← → to navigate between lessons