matrix

This article introduces a very important content in linear algebra - matrix (Matrix). It mainly explains the properties, operations and applications of matrices in homogeneous recurrence formulas with constant coefficients.

# definition

for matrix $A$, the main diagonal is $A_{i,i}$ elements.

General use $I$ To represent the identity matrix, 1 is on the main diagonal and 0 is on the other positions.

# nature

# inverse of matrix

$A$ The inverse matrix of $P$ is to make $A \times P = I$ matrix.

The inverse matrix can be found using Gaussian elimination.

# Operation

Addition and subtraction of matrices are performed element by element.

# Matrix multiplication

Matrix multiplication only makes sense if the first matrix has the same number of columns as the second matrix has the same number of rows.

Set $A$ for $P \times M$ The matrix of$B$ for $M \times Q$ The matrix of , let the matrix $C$ is a matrix $A$ with $B$ The product of

where matrix $C$ No. 1 in $i$ line of work $j$ Column elements can be expressed as:

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

If you don’t understand the above formula, it’s okay. In layman's terms, in matrix multiplication, the result $C$ The first of the matrix $i$ line of work $j$ The number of columns is determined by the matrix $A$ No. $i$ OK $M$ Numbers and matrices $B$ No. $j$ Column $M$ The numbers are obtained by multiplying them separately and then adding them together.

Matrix multiplication satisfies the associative law but does not satisfy the general commutative law.

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

In competitions, since linear recursion can be expressed in the form of matrix multiplication, matrix fast exponentiation is usually used to find an item of the linear recursion sequence.

# optimization

First, for relatively small matrices, you can consider directly manually unrolling the loop to reduce the constant.

The loop can be rearranged to improve spatial locality. Such optimizations will not change the time complexity of matrix multiplication, but will result in a constant-level improvement.

// Take the reference code below as an example
inline 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;
} // Not as good as
inline 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;
}

# Reference code

Generally speaking, a matrix can be simulated with a two-dimensional array.

struct mat { LL a[sz][sz]; inline mat() { memset(a, 0, sizeof a); } inline 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; } inline 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; } inline 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; } inline 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; }
};

# application

# Matrix accelerated recursion

Everyone should be very familiar with the Fibonacci Sequence. Among the Fibonacci numbers,$F_1 = F_2 = 1$，$F_i = F_{i - 1} + F_{i - 2}(i \geq 3)$。

If there is a question asking you to find the Fibonacci sequence number $n$ The simplest way to determine the value of an item is to directly recurse it. But if $n$ The range has reached $10^{18}$ level, recursion will not work, and TLE will be stable. Consider matrices to speed up recursion.

Set $Fib(n)$ represents a $1 \times 2$ matrix $\left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]$. We hope based on $Fib(n-1)=\left[ \begin{array}{ccc}F_{n-1} & F_{n-2} \end{array}\right]$ roll out $Fib(n)$。

Try to derive a matrix $\text{base}$,make $Fib(n-1) \times \text{base} = Fib(n)$,Right now $\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \text{base} = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]$。

How to push it? because $F_n=F_{n-1}+F_{n-2}$,so $\text{base}$ The first column of the matrix should be $\left[\begin{array}{ccc} 1 \\ 1 \end{array}\right]$, so that when performing matrix multiplication operations, $F_{n-1}$ with $F_{n-2}$ Added together, we get $F_n$. In the same way, in order to conclude $F_{n-1}$,matrix $\text{base}$ The second column of should be $\left[\begin{array}{ccc} 1 \\ 0 \end{array}\right]$。

In summary:$\text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]$ original form as $\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right] = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]$

How to convert it into code?

Define initial matrix $\text{ans} = \left[\begin{array}{ccc}F_2 & F_1\end{array}\right] = \left[\begin{array}{ccc}1 & 1\end{array}\right], \text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]$. So,$F_n$ It’s equal to $\text{ans} \times \text{base}^{n-2}$ The first row and first column elements of this matrix are $\left[\begin{array}{ccc}1 & 1\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2}$ The first row and first column elements.

Note that matrix multiplication does not satisfy the commutative law, so it must not be written as $\left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2} \times \left[\begin{array}{ccc}1 & 1\end{array}\right]$ The first row and first column elements. In addition, for $n \leq 2$ In the case of, output directly $1$ That's it, no need to perform matrix fast exponentiation.

Why multiply $\text{base}$ matrix $n-2$ power instead of $n$ What about the second power? because $F_1, F_2$ It can be found without matrix multiplication. In other words, if you only perform one multiplication, you have already found $F_3$ . If you still don’t quite understand why power is $n-2$, it is recommended to calculate by hand.

The following is the Fibonacci sequence $n$ item pair $10^9+7$ Sample code for taking modulo (core part).

const 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 complex example.

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

We found,$f_n$ and $f_{n-1}, f_{n-2}, n$ Related, so consider constructing a matrix to describe the state.

But found that if the matrix only has these three elements $\begin{bmatrix}f_n& f_{n-1}& n\end{bmatrix}$ It is difficult to construct a transfer equation because the exponentiation operation and $+1$ It cannot be described by a matrix.

So consider constructing a larger matrix.

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

We wish to construct a recursion matrix that can be transferred to

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

The transfer 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 modification

???+note ""THUSCH 2017" The Great Magician"
Produced by Little L, the great magician $n$ A magic crystal ball, each crystal ball has energy values ​​of three attributes: water, fire, and earth. Little L takes this $n$ The crystal balls are lined up on the ground from front to back, and then today's magic show begins.

We use $A_i,\ B_i,\ C_i$ to represent the energy values ​​of water, fire, and earth in the $i$-th crystal ball from front to back (the subscript starts from $1$) respectively. Little L plans to cast $m$ magic times. Each time, he will choose an interval $[l, r]$, and then cast one of the following $3$ categories and $7$ magic: 1. Magical stimulation: Make the energy of **specific attributes** in each crystal ball in the range explode, thereby enhancing the energy of another **specific attribute**. Specifically, there are three possible manifestations: - Fire element excites water element energy: Let $A_i = A_i + B_i$. - The earth element excites the fire element energy: Let $B_i = B_i + C_i$. - Water element excites earth element energy: Let $C_i = C_i + A_i$. **It should be noted that enhancing the energy of one attribute will not change the energy of another attribute. For example, $A_i = A_i + B_i$ will not increase or decrease $B_i$. ** 2. Magic enhancement: Little L waves the staff and consumes his own $v$ points of mana to change the energy of the **specific attributes** of each crystal ball in the range. Specifically, there are three possible manifestations: - Fire element energy fixed value enhancement: Let $A_i = A_i + v$. - The energy of water element is doubled and enhanced: let $B_i=B_i \cdot v$. - Earth element energy absorption and fusion: Let $C_i = v$. 3. Magic release: Little L gathers the energy of all the crystal balls in the area, fuses it into a new crystal ball, and then gives it to the audience outside the venue. The energy value of each attribute of the generated crystal ball is equal to the algebraic sum of the corresponding energy values ​​of all crystal balls in the interval. **It should be noted that the process of magic release will not actually change the energy of the crystal ball in the interval**. It is worth mentioning that the raw materials of the crystal balls manufactured and fused by Little L are all customized OI factory crystals, so these crystal balls have an energy threshold of $998244353$. When the energy value of a certain attribute in the crystal ball is greater than or equal to this threshold, the energy value will automatically modulo the threshold to prevent the crystal ball from exploding. Little W, being Little L's (only) audience member, watched the entire performance and received each of the crystal balls that Little L fused during the performance. Little W wants to know what the energy values ​​of the three attributes contained in these crystal balls are.

Since the associative and distributive laws of matrices hold, single-point modifications can be naturally extended to intervals. That is, after deriving the matrix, just use a line segment tree to maintain the interval matrix product.

A few examples will be given below.

$A_i = A_i + v$ transfer

$$\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}$$

$B_i=B_i \cdot v$ transfer

$$\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}$$

???+note ""LibreOJ 6208" Ask on the tree"
There is one tree $n$ tree of nodes, rooted at $1$ Node. Each node has two weights $k_i, t_i$, the initial values ​​are all $0$。

Three operations are given: 1. $\operatorname{Add}( x , d )$ Operation: $k_i\leftarrow k_i + d$ of all points on the path from $x$ to the root
2. $\operatorname{Mul}( x , d )$ Operation: $t_i\leftarrow t_i + d \times k_i$ of all points on the path from $x$ to the root
3. $\operatorname{Query}( x )$ operation: query the weight $t_x$ of point $x$ $n,~m \leq 100000, ~-10 \leq d \leq 10$

If you think about it directly, decentralizing operation and maintenance information is not a good idea. But matrices can be expressed 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

???+note "Problem description"
give one $n$ Order directed graph, the edge weight of each edge is $1$, and then give an integer $k$, your task is to pair all points $(u,v)$ Find out from $u$ Arrive $v$ The length is $k$ The number of paths (not necessarily simple paths, that is, points or edges on the path may be traveled multiple times).

We use this graph as an adjacency matrix $G$(For the edges in the graph $(u\to v)$,make $G[u,v]=1$, the rest are $0$ matrix; if there are multiple edges, then let $G[u,v]$ is the number of duplicate edges) represents this directed graph. The following algorithm is also applicable to graphs with self-loops.

Obviously, this adjacency matrix corresponds to $k=1$ time answer.

Suppose we know that the length is $k$ The matrix composed of the number of paths is denoted as matrix $C_k$, we want to ask $C_{k+1}$. Clearly there is a DP transfer equation

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

We can think of it as a matrix multiplication operation, so the above transfer can be described as

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

Then by expanding this recursive expression we can get

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

To calculate this matrix power, we can use the idea of ​​fast exponentiation (binary exponentiation), in $O(n^3 \log k)$ Compute the result within the complexity.

# fixed length shortest path

???+note "Problem description"
give you one $n$ Order weighted directed graph and an integer $k$. For each point pair $(u,v)$ found from $u$ Arrive $v$ of exactly contains $k$ The length of the shortest path of an edge. (It is not necessarily a simple path, that is, the points or edges on the path may be walked multiple times)

We still construct the adjacency matrix of this graph $G$，$G[i,j]$ Expresses from $i$ Arrive $j$ edge rights. if $i,j$ There is no edge between two points, then $G[i,j]=\infty$. (If there are multiple edges, the minimum value of the edge weight will be used)

Obviously the above matrix corresponds to $k=1$ The answer to the time question. we still assume we know $k$ The answer to , recorded as a matrix $L_k$. Now we want to ask $k+1$ answer. Obviously there is a transfer equation

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

In fact, we can make an analogy to matrix multiplication. You find that the above transfer only changes the sum of the products of matrix multiplication into the addition of the minimum value, so we define this operation as $\odot$,Right now

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

So get

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

Expand the recursive formula to get

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

We can still calculate the above formula using the matrix fast exponentiation method, because it is obviously associative. time complexity $O(n^3 \log k)$。

# Limited length path count/shortest path

The above algorithm is only applicable when the number of edges is fixed. However, we can improve the algorithm to solve the problem where the number of edges is less than or equal to $k$ situation. Specifically, consider the following questions:

???+note "Problem description"
give one $n$ Ordered directed graph, edge weight is $1$, and then give an integer $k$, your task is for each point pair $(u,v)$ found from $u$ Arrive $v$ length less than or equal to $k$ The number of paths (not necessarily simple paths, that is, points or edges on the path may be traveled multiple times).

Let’s simply modify this graph. We add a weight to each node as $1$ of self-loop. When walking in this way, you can walk in the loop, which is equivalent to walking in place. This includes less than or equal to $k$ situation. After modification, just do quick exponentiation of the matrix. (This algorithm still holds even if the graph has self-loops before modification).

The same method can be used to find the number of edges less than or equal to $k$ The shortest path of , that is, adding an edge weight is $0$ of self-loop.