Matrix

This article introduces a very important topic in linear algebra—matrices (Matrix). It mainly explains the properties and operations of matrices, as well as some applications of matrix multiplication.

Vectors and Matrices

In linear algebra, vectors are divided into column vectors and row vectors.

Warning

In Taiwan, China, the translations of 「column」 and 「row」 happen to be the opposite of those used in mainland China. OI Wiki follows the convention used in mainland China and uses column and row.

The main objects studied in linear algebra are column vectors, which by convention are denoted by bold lowercase letters. In linear algebra involving many vectors and matrices, when there is no risk of confusion, the vector notation above a letter may be omitted when writing by hand.

A vector is also a special kind of matrix. To represent a row vector, write a transpose symbol at the upper right of a bold lowercase letter. Row vectors generally represent equations in linear algebra.

Introduction

Matrices arise from systems of linear equations. Like vectors, matrices embody the idea of “packing” data for processing.

For example, write the system of linear equations:

$$\begin{equation} \begin{cases} 7x_1+8x_2+9x_3=13 \\ 4x_1+5x_2+6x_3=12 \\ x_1+2x_2+3x_3=11 \end{cases} \end{equation}$$

Matrices are generally represented with parentheses or square brackets. Extracting the coefficients above, write it in the form of matrix multiplication:

$$\begin{equation} \begin{pmatrix} 7 & 8 & 9 \\ 4 & 5 & 6 \\ 1 & 2 & 3 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}=\begin{pmatrix} 13 \\ 12 \\ 11 \end{pmatrix} \end{equation}$$

For short, write:

$$Ax=b$$

That is, left-multiply the column vector x of unknowns by a matrix A to obtain the column vector b. This equation can be regarded as the basic form of linear algebra.

The main operation model studied in linear algebra is the inner product. An inner product first multiplies and then adds; it is the process of left-multiplying a row vector by a column vector to obtain a number.

Matrix multiplication is an extension of the inner product. Matrix multiplication is equivalent to taking one row from the left matrix and one column from the right matrix, taking their inner product, and obtaining the corresponding element of the result matrix. The mnemonic is “left row, right column”.

When the object being studied is the column vector on the right, matrix multiplication amounts to left-multiplying the column vector. From the perspective of left multiplication, a matrix is a transformation of column vectors: transform every column vector of the matrix on the right in the matrix multiplication, obtaining the corresponding column vector in the result matrix.

A matrix can transform one column vector, or “pack” the transformation of a group of column vectors, and can even transform the whole space—that is, all column vectors. When a matrix is viewed as a transformation of the whole space, it is detached from the space and becomes a purely transformational object.

Definition

For a matrix $A$, the main diagonal refers to the elements $A_{i,i}$.

$I$ generally denotes the identity matrix, whose main diagonal consists of 1 and whose other positions are 0.

Matrices of the Same Type

Two matrices are called matrices of the same type if their numbers of rows and columns correspondingly agree.

Square Matrix

A matrix whose number of rows equals its number of columns is called a square matrix. A square matrix is a special kind of matrix. The customary phrase “an $n$-order matrix” actually refers to an $n$-order square matrix. Square matrices with the same order are matrices of the same type.

General matrices are used when studying systems of equations, vector groups, and the rank of matrices. Square matrices are used when studying eigenvalues and eigenvectors, and quadratic forms.

Main Diagonal

In a square matrix, the elements whose row and column numbers are equal form the main diagonal.

Symmetric Matrix

If the elements of a square matrix are symmetric about the main diagonal—that is, for any $i$ and $j$, the element in row $i$, column $j$ equals the element in row $j$, column $i$—then the square matrix is called a symmetric matrix.

Diagonal Matrix

A square matrix whose elements outside the main diagonal are all $0$ is called a diagonal matrix, generally written as:

$$\operatorname{diag}\{\lambda_1,\cdots,\lambda_n\}$$

The $\lambda_1,\cdots,\lambda_n$ in the expression are the elements on the main diagonal.

A diagonal matrix is a symmetric matrix.

If all elements of a diagonal matrix are $1$, it is called the identity matrix and denoted by $I$. As long as multiplication is possible, multiplying any matrix by the identity matrix leaves it unchanged, regardless of its shape.

Triangular Matrix

If all elements below the main diagonal of a square matrix are $0$, it is called an upper triangular matrix. If all elements above the main diagonal are $0$, it is called a lower triangular matrix.

The product of two upper (lower) triangular matrices is still an upper (lower) triangular matrix. If all diagonal elements are not $0$, an upper (lower) triangular matrix is invertible, and its inverse is also upper (lower) triangular.

Unit Triangular Matrix

If the diagonal of an upper triangular matrix $A$ consists entirely of $1$s, then $A$ is called a unit upper triangular matrix. If the diagonal of a lower triangular matrix $A$ consists entirely of $1$s, then $A$ is called a unit lower triangular matrix.

The product of two unit upper (lower) triangular matrices is still a unit upper (lower) triangular matrix, and the inverse of a unit upper (lower) triangular matrix is also a unit upper (lower) triangular matrix.

Operations

Linear Operations on Matrices

The linear operations on matrices are addition and subtraction, and scalar multiplication; all of them are performed element by element. Only matrices of the same type can be added or subtracted correspondingly.

Transpose of a Matrix

The transpose of a matrix is indicated by writing a transpose “T” symbol at the upper right of the matrix, meaning that the rows and columns of the matrix are exchanged.

A symmetric matrix remains unchanged before and after transposition.

Matrix Multiplication

Matrix multiplication is a generalization of the inner product of vectors.

The product of matrices is meaningful only when the number of columns in the first matrix equals the number of rows in the second matrix.

Let $A$ be a $P \times M$ matrix, $B$ be an $M \times Q$ matrix, and let matrix $C$ be the product of matrix $A$ and matrix $B$.

In matrix $C$, the element in row $i$, column $j$ can be written as:

$$C_{i,j} = \sum_{k=1}^MA_{i,k}B_{k,j}$$

In matrix multiplication, in the result matrix $C$, the number at row $i$, column $j$ is obtained from matrix $A$'s row $i$ containing $M$ numbers and matrix $B$'s column $j$ containing $M$ numbers by multiplying and then adding. Here, multiplying and then adding is the inner product of vectors. The number in row $i$, column $j$ of the product matrix is exactly the inner product of multiplier matrix $A$'s $i$-th row vector and multiplier matrix $B$'s $j$-th column vector; the mnemonic is left row, right column.

The vectors studied in linear algebra are mostly column vectors. Based on this definition of matrix multiplication, we often study left multiplication of a column vector by a matrix. It also illustrates the idea of “packed processing”: processing the inner products of many vectors at once.

Matrix multiplication is associative, but does not generally satisfy the commutative law.

Using associativity, matrix multiplication can be optimized using the idea of binary exponentiation.

In contests, since linear recurrences can be expressed in the form of matrix multiplication, matrix binary exponentiation is also commonly used to find a term of a linearly recurrent sequence.

Optimization

First, for relatively small matrices, consider manually unrolling the loops directly to reduce the constant factor.

The loops can be reordered to improve spatial locality. This optimization does not change the time complexity of matrix multiplication, but gives a constant-factor improvement.

// Taking the reference code below as an example
mat operator*(const mat& T) const {
  mat res;
  for (int i = 0; i < sz; ++i)
    for (int j = 0; j < sz; ++j)
      for (int k = 0; k < sz; ++k) {
        res.a[i][j] += mul(a[i][k], T.a[k][j]);
        res.a[i][j] %= MOD;
      }
  return res;
}

// Rather than
mat operator*(const mat& T) const {
  mat res;
  int r;
  for (int i = 0; i < sz; ++i)
    for (int k = 0; k < sz; ++k) {
      r = a[i][k];
      for (int j = 0; j < sz; ++j)
        res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
    }
  return res;
}

Inverse of a Square Matrix

For a square matrix $A$, its inverse matrix $P$ is a matrix satisfying $A \times P = I$.

An inverse matrix does not necessarily exist. If it exists, Gaussian elimination can be used to find it.

Determinant of a Square Matrix

The determinant is an operation on a square matrix.

Reference Code

Generally, a two-dimensional array can be used to simulate a matrix.

struct mat {
  LL a[sz][sz];

  mat() { memset(a, 0, sizeof a); }

  mat operator-(const mat& T) const {
    mat res;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) {
        res.a[i][j] = (a[i][j] - T.a[i][j]) % MOD;
      }
    return res;
  }

  mat operator+(const mat& T) const {
    mat res;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) {
        res.a[i][j] = (a[i][j] + T.a[i][j]) % MOD;
      }
    return res;
  }

  mat operator*(const mat& T) const {
    mat res;
    int r;
    for (int i = 0; i < sz; ++i)
      for (int k = 0; k < sz; ++k) {
        r = a[i][k];
        for (int j = 0; j < sz; ++j)
          res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
      }
    return res;
  }

  mat operator^(LL x) const {
    mat res, bas;
    for (int i = 0; i < sz; ++i) res.a[i][i] = 1;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) bas.a[i][j] = a[i][j] % MOD;
    while (x) {
      if (x & 1) res = res * bas;
      bas = bas * bas;
      x >>= 1;
    }
    return res;
  }
};

Two Perspectives on Systems of Linear Equations

There are two perspectives for viewing matrix A, or transformation A.

The first perspective is to look by rows and observe every row of A. In this way, regard A as a system of equations. This gives the process of solving equations by elimination.

The second perspective is to look by columns and observe every column of A. A itself is also composed of column vectors. This is equivalent to regarding transformation A itself as a group of column vectors, with x as the unknown coefficients, and asking whether the group of column vectors in A can be combined with the unknowns to form column vector b.

For example, the opening example becomes:

$$\begin{equation} \begin{pmatrix} 7 \\ 4 \\ 1 \end{pmatrix}x_1+\begin{pmatrix} 8 \\ 5 \\ 2 \end{pmatrix}x_2+\begin{pmatrix} 9 \\ 6 \\ 3 \end{pmatrix}x_3=\begin{pmatrix} 13 \\ 12 \\ 11 \end{pmatrix} \end{equation}$$

Solving the equations becomes the question of whether, by adjusting the three coefficients x, the given three basis vectors can be combined to form the resulting vector.

Looking by columns is more novel than looking by rows. From the column perspective, linear independence and linear dependence can be studied.

Applications of Matrix Multiplication

Matrix-Accelerated Recurrence

Take the Fibonacci sequence as an example. In the Fibonacci sequence, $F_1 = F_2 = 1$ and $F_i = F_{i - 1} + F_{i - 2}(i \geq 3)$.

If a problem asks you to find the $n$-th term of the Fibonacci sequence, the simplest method is direct recurrence. But if the range of $n$ reaches the order of $10^{18}$, recurrence is no longer feasible, so we can consider matrix-accelerated recurrence.

According to the matrix form of the recurrence formula for the Fibonacci sequence:

$$\begin{bmatrix} F_{n-1} & F_{n-2} \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} F_n & F_{n-1} \end{bmatrix}$$

Define the initial matrix $\text{ans} = \begin{bmatrix}F_2 & F_1\end{bmatrix} = \begin{bmatrix}1 & 1\end{bmatrix}, \text{base} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$. Then $F_n$ is the element in row 1, column 1 of the matrix $\text{ans} \text{base}^{n-2}$, namely the element in row 1, column 1 of $\begin{bmatrix}1 & 1\end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2}$.

Note

Matrix multiplication is not commutative, so you must not write the element in row 1, column 1 of $\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2} \begin{bmatrix}1 & 1\end{bmatrix}$. In addition, for $n \leq 2$, simply output $1$; there is no need to perform matrix binary exponentiation.

Why multiply by the $\text{base}$ matrix's $n-2$-th power rather than the $n$-th power? Because $F_1, F_2$ can be obtained without matrix multiplication. In other words, one multiplication is enough to obtain $F_3$. If it is still unclear why the exponent is $n-2$, try calculating it by hand.

Here is sample code for finding the $n$-th Fibonacci number modulo $10^9+7$ (core part).

constexpr int mod = 1000000007;

struct Matrix {
  int a[3][3];

  Matrix() { memset(a, 0, sizeof a); }

  Matrix operator*(const Matrix &b) const {
    Matrix res;
    for (int i = 1; i <= 2; ++i)
      for (int j = 1; j <= 2; ++j)
        for (int k = 1; k <= 2; ++k)
          res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
    return res;
  }
} ans, base;

void init() {
  base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
  ans.a[1][1] = ans.a[1][2] = 1;
}

void qpow(int b) {
  while (b) {
    if (b & 1) ans = ans * base;
    base = base * base;
    b >>= 1;
  }
}

int main() {
  int n = read();
  if (n <= 2) return puts("1"), 0;
  init();
  qpow(n - 2);
  println(ans.a[1][1] % mod);
}

This is a slightly more complicated example.

$$\begin{gathered} f_{1} = f_{2} = 0\\ f_{n} = 7f_{n-1}+6f_{n-2}+5n+4\times 3^n \end{gathered}$$

We see that $f_n$ is related to $f_{n-1}, f_{n-2}, n$, so we consider constructing a matrix to describe the state.

However, if the matrix contains only the three elements $\begin{bmatrix}f_n& f_{n-1}& n\end{bmatrix}$, it is difficult to construct the transition equation, because exponentiation and $+1$ cannot be described by a matrix.

Therefore, consider constructing a larger matrix.

$$\begin{bmatrix}f_n& f_{n-1}& n& 3^n & 1\end{bmatrix}$$

We want to construct a recurrence matrix that can transition to

$$\begin{bmatrix} f_{n+1}& f_{n}& n+1& 3^{n+1} & 1 \end{bmatrix}$$

The transition matrix is

$$\begin{bmatrix} 7 & 1 & 0 & 0 & 0\\ 6 & 0 & 0 & 0 & 0\\ 5 & 0 & 1 & 0 & 0\\ 12 & 0 & 0 & 3 & 0\\ 5 & 0 & 1 & 0 & 1 \end{bmatrix}$$

Matrix-Expression Updates

「THUSCH 2017」The Great Magician

The great magician Little L makes $n$ magic crystal balls, each with energy values for the water, fire, and earth attributes. Little L arranges these $n$ crystal balls in a line on the ground from front to back, and begins today's magic performance.

Let $A_i,\ B_i,\ C_i$ denote, respectively, the energy values of water, fire, and earth in the $i$-th crystal ball from the front (indices start at $1$).

Little L plans to cast $m$ spells. Each time, he chooses an interval $[l, r]$ and casts one of the following $3$ major classes and $7$ kinds of spells:

1. Magical stimulation: the energy of a specific attribute in every crystal ball in the interval erupts, thereby enhancing the energy of another specific attribute. Specifically, there are three possible forms:

   • Fire-element stimulation of water-element energy: set $A_i = A_i + B_i$.

   • Earth-element stimulation of fire-element energy: set $B_i = B_i + C_i$.

   • Water-element stimulation of earth-element energy: set $C_i = C_i + A_i$.

     Note that enhancing the energy of one attribute does not change the energy of the other attribute; for example, $A_i = A_i + B_i$ does not increase or decrease $B_i$.

2. Magical enhancement: Little L waves his staff and consumes $v$ points of his own mana to change the energy of a specific attribute in every crystal ball in the interval. Specifically, there are three possible forms:

   • Fixed-value enhancement of fire-element energy: set $A_i = A_i + v$.
   • Doubling enhancement of water-element energy: set $B_i=B_i \cdot v$.
   • Absorption and fusion of earth-element energy: set $C_i = v$.

3. Magical release: Little L gathers the energy of all crystal balls in the interval, fuses it into a new crystal ball, and gives it to an audience member outside the performance area. The energy value of each attribute in the generated crystal ball equals the algebraic sum of the corresponding energy values of all crystal balls in the interval. Note that the magical release process does not actually change the energy of the crystal balls in the interval.

It is worth mentioning that the raw materials of the crystal balls made and fused by Little L are customized OI factory crystals, so these crystal balls have an energy threshold of $998244353$. When the energy value of an attribute in a crystal ball is greater than or equal to this threshold, the energy value is automatically reduced modulo the threshold to prevent the crystal ball from exploding.

Little W is Little L's (only) audience member. He watched the entire performance and received every crystal ball fused by Little L during the performance. Little W wants to know the energy values of the three attributes contained in these crystal balls.

Because matrix multiplication satisfies associativity and distributivity, a single-point update can naturally be generalized to an interval: after deriving the matrix, maintain the product of interval matrices with a segment tree.

Here are a few examples.

The transition for $A_i = A_i + v$

$$\begin{bmatrix} A & B & C & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ v & 0 & 0 & 1\\ \end{bmatrix}= \begin{bmatrix} A+v & B & C & 1\\ \end{bmatrix}$$

The transition for $B_i=B_i \cdot v$

$$\begin{bmatrix} A & B & C & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & v & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ \end{bmatrix}= \begin{bmatrix} A & B \cdot v & C & 1\\ \end{bmatrix}$$

「LibreOJ 6208」Queries on a Tree

There is a tree with $n$ nodes, rooted at node $1$. Each node has two weights $k_i, t_i$, both initially $0$.

There are three operations:

1. $\operatorname{Add}( x , d )$ operation: for every point on the path from $x$ to the root, set $k_i\leftarrow k_i + d$.

2. $\operatorname{Mul}( x , d )$ operation: for every point on the path from $x$ to the root, set $t_i\leftarrow t_i + d \times k_i$.

3. $\operatorname{Query}( x )$ operation: query point $x$; its weight is $t_x$.

   $n,~m \leq 100000, ~-10 \leq d \leq 10$

If considered directly, it is not easy to reason about pushing operations down and maintaining information. But matrices can express this easily.

$$\begin{aligned} \begin{bmatrix}k & t & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ d & 0 & 1 \end{bmatrix} &= \begin{bmatrix}k+d & t & 1 \end{bmatrix}\\ \begin{bmatrix}k & t & 1 \end{bmatrix} \begin{bmatrix} 1 & d & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} &= \begin{bmatrix}k & t+d \times k & 1 \end{bmatrix} \end{aligned}$$

Fixed-Length Path Statistics

Problem Statement

Given an $n$-order directed graph whose every edge has weight $1$, and an integer $k$, your task is to find, for every pair of vertices $(u,v)$, the number of paths from $u$ to $v$ of length $k$ (the path need not be simple; vertices or edges on it may be traversed multiple times).

Represent this graph with its adjacency matrix $G$ (for an edge $(u\to v)$ in the graph, set $G[u,v]=1$ and all other entries to $0$; if there are multiple edges, let $G[u,v]$ be their number). The algorithm below also applies when the graph has self-loops.

Obviously, this adjacency matrix gives the answer for $k=1$.

Suppose we know the matrix formed by the numbers of paths of length $k$, denoted by $C_k$, and want to find $C_{k+1}$. The DP transition is clearly

$$C_{k+1}[i,j] = \sum_{p = 1}^{n} C_k[i,p] \cdot G[p,j]$$

We can view this as matrix multiplication, so the transition above can be described as

$$C_{k+1} = C_k \cdot G$$

Expanding this recurrence gives

$$C_k = \underbrace{G \cdot G \cdots G}_{k \text{ 次}} = G^k$$

To compute this matrix power, we can use binary exponentiation, computing the result in $O(n^3 \log k)$ time.

Fixed-Length Shortest Paths

Problem Statement

Given an $n$-order weighted directed graph and an integer $k$, find, for every pair of vertices $(u,v)$, the length of the shortest path from $u$ to $v$ containing exactly $k$ edges. (The path need not be simple; vertices or edges on it may be traversed multiple times.)

Again construct the adjacency matrix $G$ of this graph, where $G[i,j]$ represents the weight of the edge from $i$ to $j$. If there is no edge between $i,j$, then $G[i,j]=\infty$. (For multiple edges, take the minimum edge weight.)

Obviously, the matrix above gives the answer for $k=1$. Again suppose we know the answer for $k$, denoted by matrix $L_k$. We now want the answer for $k+1$. The transition is clearly

$$L_{k+1}[i,j] = \min_{1\le p \le n} \left\{L_k[i,p] + G[p,j]\right\}$$

In fact, we can compare this with matrix multiplication: the transition above merely changes the sum of products in matrix multiplication into addition followed by taking the minimum. Thus define this operation as $\odot$, namely

$$A \odot B = C~~\Longleftrightarrow~~C[i,j]=\min_{1\le p \le n}\left\{A[i,p] + B[p,j]\right\}$$

Thus

$$L_{k+1} = L_k \odot G$$

Expanding the recurrence gives

$$L_k = \underbrace{G \odot \ldots \odot G}_{k\text{ 次}} = G^{\odot k}$$

We can still compute the expression with matrix binary exponentiation, because it clearly satisfies associativity. The time complexity is $O(n^3 \log k)$.

Bounded-Length Path Counting/Shortest Paths

The algorithm above applies only when the number of edges is fixed. However, we can improve it to handle the case where the number of edges is at most $k$. Specifically, consider the following problem:

Problem Statement

Given an $n$-order directed graph with edge weights $1$ and an integer $k$, your task is to find, for every pair of vertices $(u,v)$, the number of paths from $u$ to $v$ whose length is at most $k$ (the path need not be simple; vertices or edges on it may be traversed multiple times).

For every vertex $v$, create a virtual vertex $v'$ to record the answer, and add the two edges $(v,v')$ and $(v',v')$ to the graph. Then, for a pair $(u,v)$, the number of paths from $u$ to $v$ with at most $k$ edges equals the number of paths from $u$ to $v'$ with exactly $k+1$ edges. This is because for every path with $m(m \le k)$ edges $(p_0=u)\to p_1\to p_2 \to \dots \to p_{m-1} \to (p_m=v)$, there is a corresponding path with $k+1$ edges $(p_0=u)\to p_1 \to p_2 \to \dots \to p_{m-1} \to (p_m=v) \to v' \to \dots \to v'$.

To find shortest paths with at most $k$ edges, simply add a self-loop of weight $0$ to every vertex.

Exercises

• Luogu P1962 Fibonacci Sequence, namely the example above; same problem as POJ3070
• Luogu P1349 Generalized Fibonacci Sequence, the $\text{base}$ matrix needs to be changed slightly
• Luogu P1939 Template Matrix Acceleration (Sequence), the $\text{base}$ matrix becomes a $3 \times 3$ matrix; the derivation is similar to that above.

Part of this page is translated from the blog posts Кратчайшие пути фиксированной длины, количества путей фиксированной длины and its English translation Number of paths of fixed length/Shortest paths of fixed length. The Russian version is licensed under Public Domain + Leave a Link; the English version is licensed under CC-BY-SA 4.0.