Matrix Multiplication

행렬 계산의 가장 기본적인 연산인 matrix multiplication을 중심으로 연산 구조, 계산 복잡도, 행렬 구조 활용, 메모리 효율성을 정리합니다.

---

1. Matrix and Vector Operations

1) Matrix

Keyword: Matrix

Concept: 수를 행과 열 형태로 배열한 객체.

$$A\in\mathbb{R}^{m\times n}$$

이면

$$A= \begin{bmatrix} a_{11} & \cdots & a_{1n}\\ \vdots & \ddots & \vdots\\ a_{m1} & \cdots & a_{mn} \end{bmatrix}.$$

행렬의 $(i,j)$ 성분은 $a_{ij}$로 나타냅니다.

Example Code:

#include <vector>

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

Matrix A = {
    {1.0, 2.0},
    {3.0, 4.0}
};

---

2) Vector

Keyword: Vector

Concept: 행렬의 한 열로 볼 수 있는 1차원 데이터.

$$x\in\mathbb{R}^{n}$$ $$x= \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix}.$$

Example Code:

#include <vector>

using Vector = std::vector<double>;

Vector x = {1.0, 2.0, 3.0};

---

3) Transpose

Keyword: Transpose

Concept: 행과 열을 서로 교환하는 연산.

$$C=A^T$$

이면

$$c_{ij}=a_{ji}.$$

Example Code:

Matrix transpose(const Matrix& A) {
    const std::size_t m = A.size();
    const std::size_t n = A[0].size();

    Matrix T(n, std::vector<double>(m));

    for (std::size_t i = 0; i < m; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            T[j][i] = A[i][j];
        }
    }

    return T;
}

---

2. Basic Linear Algebra Operations

1) Dot Product

Keyword: Dot Product

Concept: 두 벡터의 대응 성분을 곱한 뒤 모두 더하는 연산.

$$c=x^Ty$$ $$c=\sum_{i=1}^{n}x_i y_i.$$

Complexity:

$$O(n)$$

Example Code:

double dot(const Vector& x, const Vector& y) {
    double result = 0.0;

    for (std::size_t i = 0; i < x.size(); ++i) {
        result += x[i] * y[i];
    }

    return result;
}

---

2) Saxpy

Keyword: SAXPY

Concept: 벡터에 다른 벡터의 scalar multiple을 더하는 연산.

$$y\leftarrow ax+y.$$

SAXPY는

Scalar A X Plus Y

를 의미합니다.

많은 수치 선형대수 알고리즘을 구성하는 기본 연산입니다.

Complexity:

$$O(n)$$

Example Code:

void saxpy(double a, const Vector& x, Vector& y) {
    for (std::size_t i = 0; i < x.size(); ++i) {
        y[i] += a * x[i];
    }
}

---

3) Gaxpy

Keyword: GAXPY

Concept: matrix-vector multiplication을 포함하는 vector update.

$$y\leftarrow Ax+y.$$

각 성분은

$$y_i \leftarrow y_i+\sum_{j=1}^{n}a_{ij}x_j$$

로 계산합니다.

Complexity:

$$O(mn)$$

Example Code:

void gaxpy(
    const Matrix& A,
    const Vector& x,
    Vector& y
) {
    const std::size_t m = A.size();
    const std::size_t n = A[0].size();

    for (std::size_t i = 0; i < m; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            y[i] += A[i][j] * x[j];
        }
    }
}

---

3. Matrix-Vector Multiplication

1) Row-Oriented Multiplication

Keyword: Row-Oriented Gaxpy

Concept: 행렬의 각 행과 벡터의 dot product를 계산하는 방식.

$$y_i=A(i,:)x.$$

즉 각 행마다 하나의 dot product를 수행합니다.

Example Code:

Vector matvec_row(
    const Matrix& A,
    const Vector& x
) {
    Vector y(A.size(), 0.0);

    for (std::size_t i = 0; i < A.size(); ++i) {
        for (std::size_t j = 0; j < x.size(); ++j) {
            y[i] += A[i][j] * x[j];
        }
    }

    return y;
}

---

2) Column-Oriented Multiplication

Keyword: Column-Oriented Gaxpy

Concept: $Ax$를 행렬의 열 벡터들의 선형결합으로 해석하는 방식.

행렬을

$$A= \begin{bmatrix} a_1 & a_2 & \cdots & a_n \end{bmatrix}$$

라고 하면

$$Ax = x_1a_1+x_2a_2+\cdots+x_na_n.$$

즉 반복적인 SAXPY 연산으로 구현할 수 있습니다.

Example Code:

Vector matvec_column(
    const Matrix& A,
    const Vector& x
) {
    const std::size_t m = A.size();
    const std::size_t n = x.size();

    Vector y(m, 0.0);

    for (std::size_t j = 0; j < n; ++j) {
        for (std::size_t i = 0; i < m; ++i) {
            y[i] += A[i][j] * x[j];
        }
    }

    return y;
}

---

4. Outer Product

1) Outer Product

Keyword: Outer Product

Concept: 열 벡터와 행 벡터를 곱하여 행렬을 만드는 연산.

$$A=xy^T.$$

각 원소는

$$a_{ij}=x_i y_j.$$

Complexity:

$$O(mn)$$

Example Code:

Matrix outer(
    const Vector& x,
    const Vector& y
) {
    Matrix A(
        x.size(),
        std::vector<double>(y.size())
    );

    for (std::size_t i = 0; i < x.size(); ++i) {
        for (std::size_t j = 0; j < y.size(); ++j) {
            A[i][j] = x[i] * y[j];
        }
    }

    return A;
}

---

5. Matrix-Matrix Multiplication

1) Matrix Multiplication

Keyword: Matrix Multiplication

Concept: 행렬 $A$의 행과 행렬 $B$의 열 사이의 dot product를 계산하는 연산.

$$C=AB$$

이면

$$c_{ij} = \sum_{k=1}^{r} a_{ik}b_{kj}.$$

행렬 크기가

$$A\in\mathbb{R}^{m\times r}, \qquad B\in\mathbb{R}^{r\times n}$$

이면

$$C\in\mathbb{R}^{m\times n}.$$

Complexity:

$$O(mnr)$$

정방행렬에서는

$$O(n^3).$$

Example Code:

Matrix multiply(
    const Matrix& A,
    const Matrix& B
) {
    const std::size_t m = A.size();
    const std::size_t r = A[0].size();
    const std::size_t n = B[0].size();

    Matrix C(m, std::vector<double>(n, 0.0));

    for (std::size_t i = 0; i < m; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            for (std::size_t k = 0; k < r; ++k) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return C;
}

---

2) Dot Product Form

Keyword: Dot Product Form

Concept: 행렬 곱의 각 원소를 하나의 dot product로 계산하는 관점.

$$c_{ij} = A(i,:)B(:,j).$$

---

3) Saxpy Form

Keyword: Saxpy Form

Concept: 결과 행렬의 각 열을 $A$의 열들의 선형결합으로 계산하는 관점.

$$C(:,j) = \sum_{k=1}^{r} A(:,k)b_{kj}.$$

Example Code:

Matrix multiply_saxpy(
    const Matrix& A,
    const Matrix& B
) {
    const std::size_t m = A.size();
    const std::size_t r = A[0].size();
    const std::size_t n = B[0].size();

    Matrix C(m, std::vector<double>(n, 0.0));

    for (std::size_t j = 0; j < n; ++j) {
        for (std::size_t k = 0; k < r; ++k) {
            for (std::size_t i = 0; i < m; ++i) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return C;
}

---

4) Outer Product Form

Keyword: Outer Product Form

Concept: 행렬 곱을 여러 outer product의 합으로 표현하는 방식.

$A$의 $k$번째 열을 $a_k$, $B$의 $k$번째 행을 $b_k^T$라고 하면

$$AB = \sum_{k=1}^{r} a_kb_k^T.$$

Meaning: Matrix multiplication은 rank-1 matrix들의 합으로 볼 수 있습니다.

---

6. Loop Ordering

1) Loop Order

Keyword: Loop Ordering

Concept: Matrix multiplication의 세 반복문 $i$, $j$, $k$의 실행 순서.

가능한 순서는

• ijk
• jik
• ikj
• jki
• kij
• kji

입니다.

모두 같은 결과와 같은 수준의 flop 수를 가지지만 메모리 접근 패턴은 다를 수 있습니다.

Example Code:

ijk 방식:

for (std::size_t i = 0; i < n; ++i) {
    for (std::size_t j = 0; j < n; ++j) {
        for (std::size_t k = 0; k < n; ++k) {
            C[i][j] += A[i][k] * B[k][j];
        }
    }
}

ikj 방식:

for (std::size_t i = 0; i < n; ++i) {
    for (std::size_t k = 0; k < n; ++k) {
        for (std::size_t j = 0; j < n; ++j) {
            C[i][j] += A[i][k] * B[k][j];
        }
    }
}

Meaning: 실제 성능에서는 arithmetic count뿐 아니라 memory access pattern도 중요합니다.

---

7. Flop

1) Flop

Keyword: FLOP

Concept: 하나의 floating-point addition, subtraction, multiplication 또는 division.

대표적인 계산량은 다음과 같습니다.

Operation	Flops
$x^Ty$	approximately $2n$
$y\leftarrow ax+y$	approximately $2n$
$y\leftarrow Ax+y$	approximately $2mn$
$A\leftarrow A+yx^T$	approximately $2mn$
$C\leftarrow C+AB$	approximately $2mnr$

---

8. Big-O Complexity

1) Big-O

Keyword: Big-O

Concept: 문제 크기가 증가할 때 계산량이 얼마나 빠르게 증가하는지 나타냅니다.

주요 행렬 연산은

$$\text{Dot Product}=O(n),$$ $$\text{Matrix-Vector}=O(n^2),$$ $$\text{Matrix-Matrix}=O(n^3)$$

입니다.

---

9. BLAS Levels

1) BLAS

Keyword: BLAS

Concept: Basic Linear Algebra Subprograms.

기본적인 선형대수 연산들을 표준화한 연산 계층입니다.

---

2) Level 1 BLAS

Keyword: Level 1

Concept: Vector-vector 연산.

예:

$$x^Ty$$ $$y\leftarrow ax+y.$$

Complexity:

$$O(n)$$

---

3) Level 2 BLAS

Keyword: Level 2

Concept: Matrix-vector 연산.

예:

$$y\leftarrow Ax+y.$$

Complexity:

$$O(n^2)$$

---

4) Level 3 BLAS

Keyword: Level 3

Concept: Matrix-matrix 연산.

예:

$$C\leftarrow C+AB.$$

Complexity:

$$O(n^3)$$

Level 3 연산은 같은 데이터를 반복적으로 사용할 수 있기 때문에 고성능 matrix computation에서 특히 중요합니다.

---

10. Structured Matrices

1) Band Matrix

Keyword: Band Matrix

Concept: 주대각선 주변의 제한된 영역에만 nonzero 원소가 존재하는 행렬.

lower bandwidth가 $p$, upper bandwidth가 $q$라면

$$a_{ij}=0$$

for

$$i>j+p$$

또는

$$j>i+q.$$

Meaning: Zero 영역을 계산하거나 저장하지 않음으로써 계산량과 메모리를 줄일 수 있습니다.

---

2) Diagonal Matrix

Keyword: Diagonal Matrix

Concept: 주대각선 이외의 모든 원소가 0인 행렬.

$$D= \operatorname{diag}(d_1,\ldots,d_n).$$

Matrix-vector multiplication은 단순한 element-wise multiplication이 됩니다.

$$Dx = d\odot x.$$

Example Code:

Vector diagonal_multiply(
    const Vector& d,
    const Vector& x
) {
    Vector y(x.size());

    for (std::size_t i = 0; i < x.size(); ++i) {
        y[i] = d[i] * x[i];
    }

    return y;
}

---

3) Triangular Matrix

Keyword: Triangular Matrix

Concept: 주대각선 위 또는 아래의 성분이 모두 0인 행렬.

Upper triangular matrix:

$$a_{ij}=0 \qquad (i>j).$$

행렬 구조를 이용하면 불필요한 zero 연산을 제거할 수 있습니다.

Example Code:

Matrix multiply_upper(
    const Matrix& A,
    const Matrix& B
) {
    const std::size_t n = A.size();

    Matrix C(n, std::vector<double>(n, 0.0));

    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = i; j < n; ++j) {
            for (std::size_t k = i; k <= j; ++k) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return C;
}

---

4) Symmetric Matrix

Keyword: Symmetric Matrix

Concept: Transpose가 자기 자신과 같은 실수 정방행렬.

$$A^T=A.$$

즉

$$a_{ij}=a_{ji}.$$

따라서 전체 행렬을 저장하지 않고 삼각 영역 하나만 저장할 수 있습니다.

---

5) Permutation Matrix

Keyword: Permutation Matrix

Concept: Identity matrix의 행 또는 열 순서를 재배열한 행렬.

$$P^{-1}=P^T.$$

Permutation matrix와의 곱은 실제로 데이터 순서를 재배열하는 것과 같습니다.

Example Code:

Vector permute(
    const Vector& x,
    const std::vector<std::size_t>& p
) {
    Vector y(x.size());

    for (std::size_t i = 0; i < x.size(); ++i) {
        y[i] = x[p[i]];
    }

    return y;
}

---

11. Block Matrix

1) Block Matrix

Keyword: Block Matrix

Concept: 큰 행렬을 여러 개의 작은 submatrix로 나누어 표현하는 방식.

$$A= \begin{bmatrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{bmatrix}.$$

각 $A_{ij}$ 자체가 행렬입니다.

Meaning: Scalar 단위가 아니라 matrix block 단위로 알고리즘을 설계할 수 있습니다.

---

2) Block Matrix Multiplication

Keyword: Block Matrix Multiplication

Concept: 일반 matrix multiplication과 동일한 규칙을 block 단위로 적용합니다.

$$C_{\alpha\beta} = \sum_{\gamma} A_{\alpha\gamma}B_{\gamma\beta}.$$

Meaning: Block multiplication은 Level 3 BLAS를 많이 사용할 수 있어 고성능 계산에 유리합니다.

---

3) Blocking

Keyword: Blocking

Concept: 큰 행렬을 cache에 들어갈 정도의 작은 block으로 나누어 계산하는 방법.

Example Code:

void blocked_multiply(
    const Matrix& A,
    const Matrix& B,
    Matrix& C,
    std::size_t block_size
) {
    const std::size_t n = A.size();

    for (std::size_t ii = 0; ii < n; ii += block_size) {
        for (std::size_t kk = 0; kk < n; kk += block_size) {
            for (std::size_t jj = 0; jj < n; jj += block_size) {

                const std::size_t i_end =
                    std::min(ii + block_size, n);

                const std::size_t k_end =
                    std::min(kk + block_size, n);

                const std::size_t j_end =
                    std::min(jj + block_size, n);

                for (std::size_t i = ii; i < i_end; ++i) {
                    for (std::size_t k = kk; k < k_end; ++k) {
                        for (std::size_t j = jj; j < j_end; ++j) {
                            C[i][j] += A[i][k] * B[k][j];
                        }
                    }
                }
            }
        }
    }
}

---

12. Kronecker Product

1) Kronecker Product

Keyword: Kronecker Product

Concept: 한 행렬의 각 원소에 다른 행렬 전체를 곱하여 block matrix를 만드는 연산.

$$B\otimes C.$$

예를 들어

$$B= \begin{bmatrix} b_{11} & b_{12}\\ b_{21} & b_{22} \end{bmatrix}$$

이면

$$B\otimes C = \begin{bmatrix} b_{11}C & b_{12}C\\ b_{21}C & b_{22}C \end{bmatrix}.$$

주요 성질:

$$(B\otimes C)^T = B^T\otimes C^T,$$ $$(B\otimes C)(D\otimes F) = BD\otimes CF.$$

---

2) Vec Operation

Keyword: Vec

Concept: Matrix의 column들을 하나의 vector로 쌓는 연산.

$$\operatorname{vec}(X) = \begin{bmatrix} X(:,1)\\ X(:,2)\\ \vdots\\ X(:,n) \end{bmatrix}.$$

중요한 관계는

$$Y=CXB^T$$

와

$$\operatorname{vec}(Y) = (B\otimes C)\operatorname{vec}(X)$$

가 동등하다는 것입니다.

Meaning: Kronecker product를 직접 생성하지 않고 더 작은 matrix multiplication으로 계산할 수 있습니다.

---

13. Strassen Matrix Multiplication

1) Strassen Algorithm

Keyword: Strassen Algorithm

Concept: $2\times2$ block matrix multiplication에서 일반적인 8번의 block multiplication 대신 7번의 multiplication을 사용하는 알고리즘.

일반 matrix multiplication:

$$O(n^3)$$

Strassen multiplication:

$$O(n^{\log_2 7})$$

즉

$$O(n^{2.807\ldots}).$$

Meaning: 추가적인 덧셈을 사용하여 multiplication 횟수를 줄이는 divide-and-conquer 알고리즘입니다.

---

14. Fast Matrix-Vector Products

1) Discrete Fourier Transform

Keyword: DFT

Concept: Fourier matrix와 vector의 곱.

$$y=F_nx.$$

일반적인 matrix-vector 방식으로 계산하면

$$O(n^2)$$

입니다.

---

2) Fast Fourier Transform

Keyword: FFT

Concept: DFT matrix의 재귀적인 block structure를 이용하여 계산량을 줄이는 알고리즘.

$$O(n^2) \rightarrow O(n\log n).$$

Radix-2 FFT는 입력을 even index와 odd index로 분리하여 재귀적으로 계산합니다.

Example Code:

#include <complex>
#include <vector>
#include <numbers>

using Complex = std::complex<double>;

void fft(std::vector<Complex>& x) {
    const std::size_t n = x.size();

    if (n <= 1) {
        return;
    }

    std::vector<Complex> even(n / 2);
    std::vector<Complex> odd(n / 2);

    for (std::size_t i = 0; i < n / 2; ++i) {
        even[i] = x[2 * i];
        odd[i] = x[2 * i + 1];
    }

    fft(even);
    fft(odd);

    for (std::size_t k = 0; k < n / 2; ++k) {
        const double angle =
            -2.0 * std::numbers::pi *
            static_cast<double>(k) /
            static_cast<double>(n);

        Complex w = std::polar(1.0, angle);

        x[k] = even[k] + w * odd[k];
        x[k + n / 2] = even[k] - w * odd[k];
    }
}

---

15. Data Locality

1) Data Locality

Keyword: Data Locality

Concept: 필요한 데이터를 가능한 한 가까운 memory level에서 반복적으로 사용하는 성질.

일반적인 memory hierarchy는

$$\text{Disk} \rightarrow \text{Main Memory} \rightarrow \text{Cache} \rightarrow \text{CPU}$$

형태입니다.

Matrix computation에서는 flop 수만 줄이는 것만으로 충분하지 않습니다.

메모리에서 데이터를 가져오는 비용도 전체 성능을 크게 좌우합니다.

---

2) Unit Stride

Keyword: Unit Stride

Concept: Vector의 연속된 원소가 memory에서도 연속해서 저장되어 있는 상태.

연속적인 memory access는 일반적으로 비연속적인 access보다 효율적입니다.

Meaning: 같은 행렬 곱이라도 loop ordering에 따라 실제 실행 속도가 크게 달라질 수 있습니다.

---

3) Cache Blocking

Keyword: Cache Blocking

Concept: Matrix를 작은 block으로 나누어 cache에 유지하면서 여러 번 재사용하는 방법.

Meaning: 현대의 고성능 matrix multiplication에서는 단순 flop 수보다 data movement를 줄이는 것이 매우 중요합니다.

---

16. Parallel Matrix Multiplication

1) Parallel Matrix Multiplication

Keyword: Parallel Matrix Multiplication

Concept: Matrix multiplication 작업을 여러 processor에 나누어 수행하는 방법.

Block 단위로

$$C_{ij} = C_{ij} + \sum_k A_{ik}B_{kj}$$

를 여러 processor에 분배할 수 있습니다.

---

2) Load Balancing

Keyword: Load Balancing

Concept: 여러 processor에 계산량을 가능한 한 균등하게 분배하는 것.

한 processor에 일이 집중되면 다른 processor가 기다리게 되어 parallel efficiency가 낮아집니다.

---

3) Block-Cyclic Distribution

Keyword: Block-Cyclic Distribution

Concept: Matrix block들을 processor에 cyclic하게 분배하는 방식.

단순 block distribution보다 structured 또는 sparse matrix에서 workload를 균등하게 만들기 쉽습니다.

---

4) Communication Overhead

Keyword: Communication Overhead

Concept: Processor 사이에서 matrix block을 전달하기 위해 필요한 비용.

Parallel computation의 실제 성능은

$$\text{Computation} + \text{Communication} + \text{Synchronization}$$

에 의해 결정됩니다.

---

Essential Study Checklist

• Matrix and vector notation
• Transpose
• Dot product
• SAXPY
• GAXPY
• Matrix-vector multiplication
• Outer product
• Matrix-matrix multiplication
• Dot product / SAXPY / outer product forms
• Loop ordering
• FLOP
• Big-O complexity
• BLAS Level 1
• BLAS Level 2
• BLAS Level 3
• Band matrix
• Diagonal matrix
• Triangular matrix
• Symmetric matrix
• Permutation matrix
• Block matrix
• Blocking
• Kronecker product
• Vec operation
• Strassen algorithm
• DFT
• FFT
• Data locality
• Unit stride
• Cache blocking
• Load balancing
• Block-cyclic distribution
• Communication overhead