分类目录归档:算法

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.

线性代数相关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\) 上的一个的向量,就是要求这些向量构成矩阵的零空间的大小。

如何编写缓存友好的代码——从一道卡常数的题目说起

WC2017的第二题是一道卡常三合一题。其中,第一个子任务要求对最多2e8个32位无符号整数排序,时限3s,内存2G。输入采用伪随机数,输出采用hash,从而避免输入输出的影响。

看到题目第一反应是暴力。采用标准库中的sort函数,可以很容易地写出下面的代码

#include <cstdio>
#include <algorithm>
using namespace std;

typedef unsigned uint32_t;

inline uint32_t next_integer(uint32_t &x) {
  x = x ^ (x << 13);
  x = x ^ (x >> 17);
  x = x ^ (x << 5);
  return x;
}

const int n = 200000000;
uint32_t arr[n];

void genarray() {
  uint32_t x = 0x98765432;
  for (int i=0; i<n; i++) arr[i] = next_integer(x);
}

uint32_t output(uint32_t* ptr) {
  uint32_t ret = n * 4;
  uint32_t x = 23333333;
  for (int i=0; i<n; i++) {
    ret ^= ptr[i] + x;
    next_integer(x);
  }
  return ret;
}

int main() {
  genarray();
  sort(arr, arr+n);
  printf("%x\n", output(arr));
  return 0;
}

显然,这段代码是会TLE的。在我的电脑配置下(Intel Core i5-7200U, 8G, Ubuntu 16.04 LTS, 编译参数 g++ -O2 -m32 <source>,下同)一共跑了22.4s。

众所周知,标准库中的sort的复杂度是\(O(n \log n)\)。看来,必须使用线性时间排序才行。计数排序在这边用不了,因为32位无符号整数的取值范围太大了。那么,使用基数排序如何?将每个整数拆分为前16位和后16位,然后先以后16位为key使用计数排序,再以前16位为key使用计数排序,就可以在2倍线性时间内完成整个排序。

// ...
typedef unsigned uint32_t;
typedef unsigned short uint16_t;
typedef unsigned char uint8_t;

union {
  uint32_t u32;
  uint16_t word[2];
} arr[n], arr2[n];

uint32_t cnt[2][65536];

void counting() {
  for (int i=0; i<n; i++) {
    cnt[0][arr[i].word[0]]++;
    cnt[1][arr[i].word[1]]++;
  }
  for (int j=0; j<2; j++) partial_sum(cnt[j], cnt[j]+65536, cnt[j]);
}

void sort0() {
  for (int i=n-1; i>=0; i--)
    arr2[--cnt[0][arr[i].word[0]]].u32 = arr[i].u32;
}

void sort1() {
  for (int i=n-1; i>=0; i--)
    arr[--cnt[1][arr2[i].word[1]]].u32 = arr2[i].u32;
}

int main() {
  genarray();
  counting();
  sort0(); sort1();
  printf("%x\n", output((uint32_t*)arr));
  return 0;
}

这回运行时间为5.2s,比之前快了许多,但仍然不能在时间限制内通过。使用perf工具监测cpu用时,发现sort1和sort0函数是瓶颈,时间占比分别达到了34.37%和33.91%。再监测cache miss数量,发现程序总共发生约6.38亿次cache miss,而sort1和sort0中的占比分别为44.88%和43.74%。考虑基数排序的原理:每一趟基数排序是一次计数排序,对于相同key的元素,它们是在内存中连续放置的,而key不同的元素,则需要随机访问。通常,一趟基数排序的cache hit率可以如下估算:

\[ hit\% = \min\{1,  \frac{\#\text{cache line}}{k}\} * \left( 1 – \frac{\text{size of element}}{\text{size of cache line}}\right) \]

其中,k表示基数排序中基的大小,这里为65536。Intel i5-7200U的cache信息如下:

L1 data cache: 2 x 32 KBytes, 8-way set associative, 64-byte line size
L2 cache: 2 x 256 KBytes, 4-way set associative, 64-byte line size
L3 cache: 3 MBytes, 12-way set associative, 64-byte line size

经过简单的计算发现L3 cache一共只有约5万行,也就是说,这样写连3级缓存都不能充分利用,CPU必须经常花费数倍的时间访问主存,因此效率十分低下。

L1 cache一共只有512行,那么我们将32位整数分为4段,每段8位,然后进行基数排序行不行呢?让我们重新改写一下程序。

// ...

uint32_t cnt[4][256];

void counting() {
  for (int i=0; i<n; i++) {
    cnt[0][arr[i].byte[0]]++;
    cnt[1][arr[i].byte[1]]++;
    cnt[2][arr[i].byte[2]]++;
    cnt[3][arr[i].byte[3]]++;
  }
  for (int j=0; j<4; j++)
    partial_sum(cnt[j], cnt[j]+256, cnt[j]);
}

void sort0() {
  for (int i=n-1; i>=0; i--)
    arr2[--cnt[0][arr[i].byte[0]]].u32 = arr[i].u32;
}

void sort1() {
  for (int i=n-1; i>=0; i--)
    arr[--cnt[1][arr2[i].byte[1]]].u32 = arr2[i].u32;
}

void sort2() {
  for (int i=n-1; i>=0; i--)
    arr2[--cnt[2][arr[i].byte[2]]].u32 = arr[i].u32;
}

void sort3() {
  for (int i=n-1; i>=0; i--)
    arr[--cnt[3][arr2[i].byte[3]]].u32 = arr2[i].u32;
}

int main() {
  genarray();
  counting();
  sort0(); sort1(); sort2(); sort3();
  printf("%x\n", output((uint32_t*)arr));
  return 0;
}

跑下来运行时间为2.7s,能够在时限内通过。使用性能分析工具perf,发现总共只发生了2.5亿次cache miss。虽然经过修改后,理论运行时间变为2倍,但实际运行时间反而减少了将近一半。

CPU在访问主存时,往往需要等待数百个时钟周期。由于burst机制的存在,CPU可以一次性与主存交换大量连续数据,使得cache能够被十分有效的实现。要增加cache的利用率,主要是要提高程序访问数据的空间局部性和时间局部性。将 \(O(n \log n)\) 的排序改为线性时间排序,是一种典型的space-time tradeoff, 如果排序的基数太大,那么count数组会变得很大,排序后的数据在内存中的分布会很散,降低了缓存利用率,这时就可以通过适当降低基数的方式,提高cache利用率,从而达到理论运行时间和cache利用率之间的平衡。

k次幂前缀和的线性做法

给定正整数 \(x, k\) 和质数 \(p\),求
\[ \sum_{i=1}^{x} i^k \bmod p \]

这是一道很经典的题目,在杜瑜皓的《多项式及求和》中有详细讨论。比较暴力的做法是对于每个 \(i^k\) 都用快速幂算出其值,然后相加,时间复杂度为 \(O(x \log k) \),当 \(x\) 很大时,就超时了。然而利用拉格朗日插值法,这题最快可以做到 \(O(k)\) 的复杂度,其做法主要由以下几部分构成。

线性预处理逆元

求 1..n 之间每个数的逆元,如果都用费马小定理或者扩展欧几里得算,那么复杂度将会达到 \(O(n \log p)\) 或 \(O(n \log n)\)。利用一些递推式,可以线性地求出 1..n 中每个数的逆元,从而复杂度可以减少一个 log。常用的一个递推公式是

\[ i^{-1} = -\left\lfloor \frac{p}{i} \right\rfloor \cdot (p \bmod i)^{-1} \bmod p \]

这一公式的正确性也很好证明,等式右端可以写成

\[ -\frac{\left\lfloor \frac{p}{i} \right\rfloor}{p-\left\lfloor \frac{p}{i} \right\rfloor i} = \frac{1}{i} \]

利用此公式,可以在 \(O(n)\) 时间内求出 1..n 中每一个数的逆元。

线性筛求积性函数的值

注意到 \(i^k\) 是一个积性函数,这样,利用线性筛,我们只需要用快速幂计算素数 \(i\) 的值,其他点的函数值就可以直接得到。根据素数定理,前 \(n\) 个数中素数的个数约为 \(O(\frac{n}{\log n})\),从而利用快速幂计算出素数处的值,可以做到 \(O(\frac{n}{\log n}\log{k}\) ,再加上筛法本身 \(O(n)\) 的复杂度,我们可以在 \(O(n + \frac{n}{\log n}\log{k})\) 的时间复杂度内计算出 \(i^k\) 的前 \(n\) 项值。

拉格朗日插值法

注意到答案一定是关于 \(x\) 的 \(k+1\) 次多项式。而根据线性代数的知识,我们只要知道这个多项式的任意 \(k+2\) 个点处的值,我们就能确定这个多项式。拉格朗日插值公式就给出了这个多项式

\[ f(x) = \sum_{i=0}^{k+1} y_i l_i(x) \]

其中

\[ l_i(x) = \prod_{0 \leq m \leq k+1, m \neq i} \frac{x – x_m}{x_j – x_m} \]

具体到这道题上,插值公式就是

\[ f(n) = \sum_{i = 0}^{k+1} (-1)^{k+i+1} f(i) \frac{n(n-1) \cdots (n-i+1)(n-i-1) \cdots (n-d)}{(d-i)!i!} \]

其中 \(f(i)\) 就是 \(i^k\) 的前缀和,上式右端求和的每一项,分子都是 n 到 n-d 的前缀积乘后缀积,分母可以用递推线性求出,这样就可以在 \(O(k)\) 时间内解决整个问题。

#include <bits/stdc++.h>
using namespace std;

#define rep(i, n) for (int i=0; i<(n); i++)
#define Rep(i, n) for (int i=1; i<=(n); i++)
#define range(x) (x).begin(), (x).end()
typedef long long LL;
typedef unsigned long long ULL;

#define pow owahgrauhgso

ULL n, k, m;

ULL powmod(ULL b, ULL e) {
  ULL r = 1;
  while (e) {
    if (e & 1) r = r * b % m;
    b = b * b % m;
    e >>= 1;
  }
  return r;
}

const int MAXN = 1000006;
ULL inv[MAXN];

void init_inv() {
  inv[1] = 1;
  for (int i = 2; i <= k+1; i++) { 
    inv[i] = (m - m / i) * inv[m % i] % m; 
    assert(inv[i] * i % m == 1);
  }
}

ULL pow[MAXN];
ULL prime[MAXN], cnt;

void sieve() {
  pow[1] = 1;
  for (int i = 2; i <= k+1; i++) {
    if (!pow[i]) {
      pow[i] = powmod(i, k);
      prime[cnt++] = i;
      for (int j = 0; j < cnt && i*prime[j] <= k+1; j++) {
        pow[i * prime[j]] = pow[i] * pow[prime[j]] % m;
        if (i % prime[j] == 0) break;
      }
    }
  }
}

ULL sum[MAXN];
ULL ans[MAXN];

auto addmod = [](ULL a, ULL b) -> ULL {return (a+b)%m;};

ULL lagrange() {
  ULL p;
  p = 1;
  for (int i=0; i<=k+1; i++) {
    if (i) p = p * inv[i] % m;
    ans[i] = (k+1-i)&1 ? m-sum[i] : sum[i];
    ans[i] = ans[i] * p % m;
    p = p * (m + n - i) % m;
  }
  p = 1;
  for (int i=k+1; i>=0; i--) {
    if (k+1-i) p = p * inv[k+1-i] % m;
    ans[i] = ans[i] * p % m;
    p = p * (m + n - i) % m;
  }
  return accumulate(ans, ans+k+2, 0, addmod); 
}

 

int main() { cin >> n >> k >> m; init_inv(); sieve(); partial_sum(pow, pow+k+2, sum, addmod); cout << lagrange() << endl; return 0; }