标签归档:线性代数

Solving Sparse Linear System in O(nm) Time

There are many practical applications that involve solving sparse linear system \( \mathbf{A}\mathbf {x} = \mathbf{b} \), such as random walk in sparse graphs. Traditional techniques such as Gauss elimination or LU decomposition applies for general linear systems, and they run in cubic time. Fastest known algorithms for general matrix run in about \( O(n^{2.373}) \) time. However, for sparse matrix, we may achieve a much faster result. Specifically, we may solve the system in \(O(mn)\) time, where \(m\) is the number of nonempty entries in \(\mathbf{A}\).  We assume that the matrix \( \mathbf{A} \) is non-singular (i.e., invertible) throughout this passage.

1. Inverting a matrix from its annealing polynomial

Given polynomial \( p(x) \), if \( p(\mathbf{A}) = \mathbf{0} \) for matrix \( \mathbf{A} \), then we say \(p(x)\) is an annealing polynomial for matrix \( \mathbf{A} \). Due to Hamilton-Cayley theorem, the characteristic polynomial \(\det( x\mathbf{I} – \mathbf{A} ) \) is an annealing polynomial for \( \mathbf{A} \). Among all annealing polynomials for matrix \( \mathbf{A} \), the one with minimum degree is called the minimal polynomial for \( \mathbf{A} \). It can be proved that the minimum polynomial for given matrix is unique up to a constant factor, and any other annealing polynomial is a polynomial multiple of the minimum polynomial.

We can invert a matrix \( \mathbf{A} \) if we have a minimal polynomial for \( \mathbf{A} \) :
\[ a_0 I + a_1  \mathbf{A} + a_2 \mathbf{A}^2 + \cdots + a_k \mathbf{A}^k = 0 \tag{*} \]Since \( \mathbf{A} \) is invertible, 0 is never a root of its characteristic polynomial, hence we must have \(a_0 \neq 0 \). Multiplying \( \mathbf{A}^{-1} \) on both sides yields
\[ \mathbf{A}^{-1} = -\frac{a_1 I + a_2 \mathbf{A} + a_3 \mathbf{A}^2 + \cdots + a_k \mathbf{A}^{k-1}}{a_0} \]this means that we may represent the inversion of \( \mathbf{A} \) by the linear combination of powers of  \( \mathbf{A} \).

2. The Berlekamp-Massey algorithm

Berlek-Massey algorithm solves the following problem in $O(n^2)$ time:

Given a finite sequence \(\{x_i\}_{i=1}^n\), find a minimum order linear recurrence consistent with the given sequence. Formally, find a shortest sequence \(c_0 = 1, c_1, \cdots, c_{k-1}\), such that \(\sum_{l=0}^{k-1} x_{j-l} c_l = 0\) holds for all possible \(j\).

 This algorithm has many real world applications. The most typical one is to find the shortest linear feedback shift register for a given binary sequence. Also, it can be viewed as interpolating a sequence with exponential terms. One important fact for Berlekamp-Massey algorithm is, for an order-\(r\) linearly recurrent sequence, taking the first \(2r\) elements as the input of the algorithm suffices to recover the recurrence.

3. Finding the minimum polynomial 

Note that the annealing polynomial is exactly the linear recurrence of powers of a matrix. However, it is infeasible to compute the minimum polynomial from the powers of \(\mathbf{A}\). However, we may randomly pick vectors \(\mathbf{u}\) and \(\mathbf{v}\) and compute the minimum polynomial from \(\mathbf{u}^T \mathbf{A}^i \mathbf{v}\). We claim without proof that, with high probability, the coefficients of the recurrence of the sequence are exactly those of the minimum polynomial. The sequence can be computed in \(O(mn)\) time by iteratively doing sparse matrix-vector multiplication in $O(m)$ time. Finally, apply Berlekamp-Massey algorithm to the given sequence.

4. Solving the linear system

Since the inverse of a sparse matrix is generally not sparse, we won’t actually compute the inverse of \(\mathbf{A}\).  Actually, we can compute \(\mathbf{A}^{-1}\mathbf{b}\) via formula (*) in \(O(mn)\) time. The procedure is exactly the same as in finding minimum polynomial: just iteratively perform spare matrix-vector multiplication.

[Codeforces Round #513] H. Sophisticated Device

题目描述

定义一种机器,一共有5000个寄存器,前两个寄存器初始值分别为a, b,其他寄存器初始值均为1。该机器支持两种操作:

  1. ADD %dest, %src1, %src2:R[dest] = R[src1] + R[src2]
  2. POW %d, %s:R[dest] = pow(R[src], d)

所有操作均模p。要求用不超过5000条指令计算a*b。

题解

首先有加法我们就可以用类似快速幂的方式实现乘一个常数。这样,我们就可以构造出常数0(只需要乘p即可),实现减法(乘p-1获得相反数),除以一个常数(乘上它的逆元)。

精彩的是,我们可以用上述两条指令计算一个数的平方。考虑\(x^d, (x+1)^d, \cdots, (x+d)^d\)的二项展开,每一个都可以表示成\(1, x, x^2, \cdots, x^d\)的线性组合。这样,通过矩阵求逆就可以用\(x^d, (x+1)^d, \cdots, (x+d)^d\)表示出\(x^2\),从而完成了平方的计算。

有了平方操作,就可以用\(xy = ((x+y)^2-x^2-y^2)/2\)计算出两数之积了。

2018牛客暑期ACM多校训练营(第一场)题解

A. Monotonic Matrix

给定\(1\leq m,n \leq 1000\),求满足下列条件的\(m \times n\)矩阵个数,使得

  • \(a_{i,j} \in \{0, 1, 2\}\)
  • 每一行从左到右和每一列从上到下单调不减

题解:考虑每一行大于等于1和大于等于2的位置,其实就是要求两个单调不减的数列\(\{a_i\}, \{b_i\}\)的个数,满足对于任意\(i\),有\(a_i \leq b_i\)。

B. Symmetric Matrix

给定\(1 \leq n \leq 10^5\),求满足每个元素为非负整数,主对角线全0,每一行和为2的\(n \times n\)的对称矩阵的数量。

题解:把矩阵想象成邻接矩阵,就是要求顶点带标号的所有不允许自环、但允许重边的2-正则图的数量。

C. Fluorescent 2

在模2域下,给定\(m \times n\)的01矩阵\(M\),\(1 \leq n \leq 50\),\(1 \leq m \leq 1000\),求对于每一个\(n\)维向量\(v\),\(Mv\)中0的个数的三次方的和。

D. Two Graphs

给定图\(G_1 = (V, E_1), G_2 = (V, E_2)\),求有多少\(E_1\)的子集\(E’\),使得\((V, E’)\)与\(G_2\)同构。点数不超过8。

题解:枚举所有的排列,如果构成子图,就把选出来的边用bitmap表示出来,最后去重即可。

E. Removal

给定\(1 \leq n \leq 10^5\)个数的序列,每个数为不超过\(k\)的正整数。问删掉\(m\)个数后,有多少种不同的子序列。\(k, m\)不超过10。

题解:用dp[i][j][k]表示前i个数删去j个数后,最后一个数为k的方案数。可以在\(O(1)\)时间内转移。

F. Sum of Maximum

给定数列\(a_1, a_2, \cdots, a_n\),求

\[\sum_{x_1=1}^{a_1} \sum_{x_2=1}^{a_1} \cdots \sum_{x_n=1}^{a_n} \max\{x_1, x_2, \cdots, x_n\}\]

题解:把\(\max\{x_1, x_2, \cdots, x_n\}\)写为\(\sum_{m=1}^a [\exists i, x_i \geq m]\),进而改写为\(\sum_{m=1}^a 1-[\forall i, x_i < m]\)。将\(a_i\)排序后,从大到小分段计算,每一段都是\(x_i\)的一个前缀积乘以\(i^k\)的求和,后者可用Lagrange插值法快速计算。

G. Steiner Tree

H. Longest Path

I. Substring

给一个字符集为{a,b,c}的字符串,求在同构意义下的子串的等价类的数目。两个串同构当且仅当长度相等,并且存在字符集到字符集的双射f,使得一个串经过f变换后变成另一个串。

题解:将原串经过6种双射f后变成的6个串丢到SAM里,求出所有本质不同的子串。含有大于等于2种子串有6中不同的同构串,而只含有1种字符的串只有3中不同的同构串,分类讨论一下即可。

J. Different Integers

给一个序列,每次询问给出一个区间,求序列扣除这个区间后,有多少个不同数。

题解:在某次询问中,某数不出现当且仅当这个数第一次出现到最后一次出现包含于询问区间中即可。只要离线排序一下,按右端点的顺序处理即可。

线性代数相关OJ题选讲

这篇文章是关于线性代数的OJ题选讲,这几道题目都不难,主要是考察思维的灵活程度。


Fast Matrix Calculation

2014 Multi-University Training Contest 9, available at HDUOJ 4965

给定 \(N\times k\) 的矩阵 \(A\),\(k \times N\) 的矩阵 \(B\) \((2 \leq k \leq 6, 4 \leq N \leq 1000)\):
Step 1. 计算 \(N \times N\) 的矩阵 \(C = AB\).
Step 2. 计算 \(M = C^{N \times N} \bmod 6\).
Step 3. 输出 \(M\) 中所有元素的和.

题解

根据矩阵乘法的结合律,

\[M = (AB)^{N \times N} = A (BA)^{N \times N – 1} B\]

时间复杂度由裸的快速幂\(O(N^3 \log N)\)变成了\(O(N^2k + Nk^2 + k^3 \log N)\)。


Matrix Multiplication

2014 Multi-University Training Contest 5, available at HDUOJ 4920

给定\(n \times n\)的矩阵\(A, B\),求模3意义下的乘积。 \((1 \leq n \leq 800)\)

题解

令 \(M_x\): \(M_x(i, j) = [M(i, j) = x]\).

\[ AB = (A_1 – A_2)(B_1 – B_2) = A_1B_1 – A_1B_2 – A_2B_1 + A_2B_2 \]

其中 \(A_1, A_2, B_1, B_2\) 是 01矩阵. 使用 std::bitset,可以做做到\(O(4n^3/wordsize)\)。

(当然,使用其他一些优秀的卡常方法也能过这题)


Matrix Revolution

HDUOJ 1759

给定一个 \(n \times n\) 的矩阵 \(A\),计算下述矩阵的非零元素的数量

\[ A + A^2 + \cdots + A^k\]

\(A\) 是稀疏矩阵,用 \(m\) 个三元组 \((a, b, c)\) 表示,每个三元组代表 \(A_{ab} = c\)。此外,保证主对角线上的元素全为 1。其他元素均为0。

\((0 \leq n \leq 1000, 0 \leq m \leq 10n, n < k < 10^{100}, 0 < c < 10^9)\)

题解

利用图的邻接矩阵表示,可知答案就是距离不超过 \(k\) 的点对数。


Matrix Multiplication

POJ 3318

给定 \(n \times n\) 的矩阵 \(A, B, C\)。判断\(A B = C\)是否成立。\((n \leq 500, |A_{ij}|, |B_{ij}| \leq 100, |C_{ij}| \leq 10,000,000)\)

提示:有多组测试数据,\(O(n^3)\) 的算法会 TLE。

题解

随机生成 \(n\) 维列向量 \(\vec{x}\),判断 \(AB\vec{x} = C\vec{x}\) 是否成立。可在 \(O(n^2)\) 时间内完成。

算法出错当且仅当 \((AB-C)\vec{x} = 0\) 而 \(AB – C \neq 0\)。 假设我们在 \(\mathbb{F}_p\) 里做。令 \(M = AB – C\), 则 \(\vec{x}\) 在 \(M\) 的零空间中。仔细想一想 \(M\) 的零空间最多 \(n-1\) 维,因此出错的概率不会超过\(1/p\)。

(当然,使用三次方的算法加一些优秀的卡常,也许也能过)

Remark: 

这一算法被称为Freivalds算法,是一种单边错Monte Carlo随机算法。对于整数矩阵,存在确定性的\(O(n^2)\)的算法判断乘法是否正确。不过由于Freivalds算法非常易于实现,目前实际上仍然大量使用该算法 。


Square

UVA 11542, also available at vjudge

给定 \(n\) 个整数的集合你可以生成 \(2^n-1\) 个非空子集。求有多少个子集,其中的数之积为完全平方数。

\(1 \leq T \leq 30\), \(1 \leq n \leq 100\), 所有整数在 1 和 \(10^{15}\) 之间。整数的素因子最多为 500。

题解

每个整数的唯一分解序列都看成 \(\mathbb{F}_2\) 上的一个的向量,就是要求这些向量构成矩阵的零空间的大小。