矩阵在计算机程序中的应用
先總結一下矩陣運算的模版程序
求逆的方法目前不懂，求路過的大神解釋。
/* 這部分請忽略
矩陣相乘有分治來優化的方法，Strassen算法，矩陣可以填0的方法計算令它成為2^n * 2^n，貼一下分治的那部分
void Strassen(int n, T A[][N], T B[][N], T C[][N]) {T A11[N][N], A12[N][N], A21[N][N], A22[N][N];T B11[N][N], B12[N][N], B21[N][N], B22[N][N];T C11[N][N], C12[N][N], C21[N][N], C22[N][N];T M1[N][N], M2[N][N], M3[N][N], M4[N][N], M5[N][N], M6[N][N], M7[N][N];T AA[N][N], BB[N][N];if (n == 2) { //2-order Matrix_Multiply(A, B, C);} else {//將矩陣A和B分成階數相同的四個子矩陣，即分治思想。for (int i = 0; i < n / 2; i++) {for (int j = 0; j < n / 2; j++) {A11[i][j] = A[i][j];A12[i][j] = A[i][j + n / 2];A21[i][j] = A[i + n / 2][j];A22[i][j] = A[i + n / 2][j + n / 2];B11[i][j] = B[i][j];B12[i][j] = B[i][j + n / 2];B21[i][j] = B[i + n / 2][j];B22[i][j] = B[i + n / 2][j + n / 2];}}//Calculate M1 = (A0 + A3) × (B0 + B3)Matrix_Add(n / 2, A11, A22, AA);Matrix_Add(n / 2, B11, B22, BB);Strassen(n / 2, AA, BB, M1);//Calculate M2 = (A2 + A3) × B0Matrix_Add(n / 2, A21, A22, AA);Strassen(n / 2, AA, B11, M2);//Calculate M3 = A0 × (B1 - B3)Matrix_Sub(n / 2, B12, B22, BB);Strassen(n / 2, A11, BB, M3);//Calculate M4 = A3 × (B2 - B0)Matrix_Sub(n / 2, B21, B11, BB);Strassen(n / 2, A22, BB, M4);//Calculate M5 = (A0 + A1) × B3Matrix_Add(n / 2, A11, A12, AA);Strassen(n / 2, AA, B22, M5);//Calculate M6 = (A2 - A0) × (B0 + B1)Matrix_Sub(n / 2, A21, A11, AA);Matrix_Add(n / 2, B11, B12, BB);Strassen(n / 2, AA, BB, M6);//Calculate M7 = (A1 - A3) × (B2 + B3)Matrix_Sub(n / 2, A12, A22, AA);Matrix_Add(n / 2, B21, B22, BB);Strassen(n / 2, AA, BB, M7);//Calculate C0 = M1 + M4 - M5 + M7Matrix_Add(n / 2, M1, M4, AA);Matrix_Sub(n / 2, M7, M5, BB);Matrix_Add(n / 2, AA, BB, C11);//Calculate C1 = M3 + M5Matrix_Add(n / 2, M3, M5, C12);//Calculate C2 = M2 + M4Matrix_Add(n / 2, M2, M4, C21);//Calculate C3 = M1 - M2 + M3 + M6Matrix_Sub(n / 2, M1, M2, AA);Matrix_Add(n / 2, M3, M6, BB);Matrix_Add(n / 2, AA, BB, C22);//Set the result to C[][N]for (int i = 0; i < n / 2; i++) {for (int j = 0; j < n / 2; j++) {C[i][j] = C11[i][j];C[i][j + n / 2] = C12[i][j];C[i + n / 2][j] = C21[i][j];C[i + n / 2][j + n / 2] = C22[i][j];}}} }算法核心就是把每次相乘的8次矩陣相乘變成了7次乘法，看別小看這減少的1次乘法，因為每遞歸1次，性能就提高了1/8，比如對于1024*1024的矩陣，第1次先分解成7次512*512的矩陣相乘，對于512*512的矩陣，又可以繼續遞歸分解成256*256的矩陣相乘，…，一直遞歸下去，假設分解到64*64的矩陣大小后就不再遞歸，那么所花的時間將是分塊矩陣乘法的(7/8) * (7/8) * (7/8) * (7/8) = 0.586倍，提高了快接近一倍
*/
?
回歸正題，看O(n^3)的矩陣乘法，這個在程序里就很可以了。
#include <iostream> #include <cmath> using namespace std;#define MAXN 4 #define zero(x) (fabs(x)<1e-8)struct mat {int n, m;double data[MAXN][MAXN];mat(){memset(data, 0, sizeof(data));} }; void debug(mat tmp) {for (int i = 0; i < tmp.n; i++) {for (int j = 0; j < tmp.m; j++)printf("%.3lf ", tmp.data[i][j]);puts("");}puts(""); } mat mul(mat& a, mat& b) { //a乘b = cint i, j, k;mat c;c.n = a.n, c.m = b.m;for (i = 0; i < c.n; i++)for (j = 0; j < c.m; j++) {for (c.data[i][j] = k = 0; k < a.m; k++)c.data[i][j] += a.data[i][k] * b.data[k][j]; //c.data[i][j] = (c.data[i][j] + ((LL) a.data[i][k] * b.data[k][j]) % mod) % mod;if (zero(c.data[i][j]))c.data[i][j] = 0;}return c; } mat pow(mat a, int x) {mat b;b.m = b.n = max(a.m, a.n);for (int i = 0; i < b.n; i++)b.data[i][i] = 1;while (x) {if (x & 1) {b = mul(a, b);}a = mul(a, a);x >>= 1;}return b; } //二分計算A^l+...A^r mat cal(mat A, int l, int r) {if (l == r)return pow(A, l);int m = (l + r + 1) / 2;mat t = cal(A, l, m - 1); //算前一半的和mat t2 = pow(A, m - l); //中介的矩陣t = add(t, mul(t2, t));if ((r - l + 1) & 1)t = add(t, pow(A, r)); //如果是奇數，則加上最后一個return t; } mat inv(mat a) { // 高斯—約旦法求逆矩陣int i, j, k, is[MAXN], js[MAXN];double t;for (k = 0; k < a.n; k++) { // 從第k行k列開始右下角子陣中選取絕對值最大的元素 記住元素所在行列號for (t = 0, i = k; i < a.n; i++)for (j = k; j < a.n; j++)if (fabs(a.data[i][j]) > t)t = fabs(a.data[is[k] = i][js[k] = j]);if (is[k] != k) // 在通過行交換和列交換將它交換到主元素位置上for (j = 0; j < a.n; j++)swap(a.data[k][j], a.data[is[k]][j]);if (js[k] != k)for (i = 0; i < a.n; i++)swap(a.data[i][k], a.data[i][js[k]]);a.data[k][k] = 1 / a.data[k][k];for (j = 0; j < a.n; j++)if (j != k)a.data[k][j] *= a.data[k][k];for (i = 0; i < a.n; i++)if (i != k)for (j = 0; j < a.n; j++)if (j != k)a.data[i][j] -= a.data[i][k] * a.data[k][j];for (i = 0; i < a.n; i++)if (i != k)a.data[i][k] *= -a.data[k][k];}for (k = a.n - 1; k >= 0; k--) {for (j = 0; j < a.n; j++)if (js[k] != k)swap(a.data[k][j], a.data[js[k]][j]);for (i = 0; i < a.n; i++)if (is[k] != k)swap(a.data[i][k], a.data[i][is[k]]);}return a; }double det(mat& a) { //求行列式int i, j, k, sign = 0;double b[MAXN][MAXN], ret = 1, t;if (a.n != a.m)return 0;for (i = 0; i < a.n; i++)for (j = 0; j < a.m; j++)b[i][j] = a.data[i][j];for (i = 0; i < a.n; i++) {if (zero(b[i][i])) {for (j = i + 1; j < a.n; j++)if (!zero(b[j][i]))break;if (j == a.n)return 0;for (k = i; k < a.n; k++)t = b[i][k], b[i][k] = b[j][k], b[j][k] = t;sign++;}ret *= b[i][i];for (k = i + 1; k < a.n; k++)b[i][k] /= b[i][i];for (j = i + 1; j < a.n; j++)for (k = i + 1; k < a.n; k++)b[j][k] -= b[j][i] * b[i][k];}if (sign & 1)ret = -ret;return ret; }int main() {freopen("data3.txt", "r", stdin);int n, m;while (~scanf("%d%d", &n, &m)) {mat ma;ma.n = n, ma.m = m;for (int i = 0; i < n; i++)for (int j = 0; j < m; j++)scanf("%lf", &ma.data[i][j]);mat ni = inv(ma); // debug(ni);puts("");mat mu = mul(ma, ni); // debug(mu); }return 0; }?
1.數學應用和圖論應用（傳遞閉包）
先來說一下傳遞閉包的問題，就是離散數學中二元關系傳遞轉化為矩陣相乘
比如POJ 3660這題，通過矩陣乘得到所有的有關系的二元點對，這題要求多少個點排名能夠確定，就是知道這個點和所有其他點的勝負關系，轉化為點的入度與出度和是否為n-1的問題
這題要注意，矩陣乘循環的順序
#include <iostream> #include <cstdio> using namespace std;#define MAXN 110 int n; bool g[MAXN][MAXN]; void mul() {int i, j, k;for (k = 1; k <= n; k++)for (i = 1; i <= n; i++)for (j = 1; j <= n; j++)g[i][j] = g[i][j] | (g[i][k] & g[k][j]); } int main() { // freopen("data4.txt", "r", stdin);int m, u, v;scanf("%d%d", &n, &m);memset(g, 0, sizeof(g));while (m--) {scanf("%d%d", &u, &v);g[u][v] = 1;}mul();int cnt = 0;for (int i = 1; i <= n; i++) {int cnt_t = 0;for (int j = 1; j <= n; j++)if (g[i][j] || g[j][i]) {cnt_t++;}if (cnt_t == n - 1)cnt++;}cout << cnt << endl;return 0; }?
1?給定一個有向圖，問從A點恰好走k步（允許重復經過邊）到達B點的方案數mod p的值
2 給出一張無向連通圖，求S到E經過k條邊的最短路
?
2.通過構造解決數學問題
我們先來考慮Fibonacci數列的遞推公式，f[i] = f[i-1] + f[i-2]
我們可以用一個轉移矩陣來表示這種轉移
有一個矩陣 ?1 1 ?乘上 ?f[i-1] 得到 f[i]
?1 0 f[i-2] ?f[i-1]
這樣f[i-1]通過一個轉移矩陣遞推出f[i]
同理如 f[i] = a1*f[i-1] + a2*f[i-2] + ~~ + ak*f[i-k]
就可以找到矩陣 a1 a2 a3 ~~ ak 通過第一行來遞增，通過之后的行來保留原來的值
1 ? 0 ?0 ? ~~ ?0
0 ? 1 ?0 ? ~~ ?0
| ? ?| ? | ? ~~ ?|
0 ? 0 ~~ ? 1 ? 0
則 f[n] = 轉移矩陣^n * f[0] 更高效的求出 f[n]
另一個應用 矩陣表示向量的縮放，翻轉，旋轉，平移。
縮放，翻轉，旋轉 可以給坐標矩陣 x 乘上一個2*2矩陣得到，如縮放乘 k ?0 翻轉乘 1 ?0 ?旋轉乘 cosa ?-sina
?y ? ? 0 ?k ?0 ?-1 ?sina cosa
但是不能滿足平移操作的要求
我們的解決辦法是加一個維度，用 x 來表示坐標，平移操作乘上矩陣 1 ?0 ?x
y ?0 ?1 ?y
1 ?0 ?0 ?1
保證統一性，其他操作都變為3*3的矩陣，3行3列為1，其他補0。
找到一個很好地圖說明
?
另外也可以解決對序列的置換
?
構造矩陣例題 HDU 1757
f(x) = a0 * f(x-1) + a1 * f(x-2) + a2 * f(x-3) + …… + a9 * f(x-10);
|f(10) | ? ? |a0 a1 a2 ...a8 a9| |f(9)|
| f(9) ?| ? ? | 1 0 0 ... 0 0| ? ? ? ?|f(8)|
| ?.... ? | = ?|.. ?... ... ?... ..| ? ? ? | .. |
| f(2) ?| ? ? | 0 0 0 ... 0 0| ? ? ? ?|f(1)|
| f(1) ?| ? ? | 0 0 0 ... 1 0| ? ? ? ?|f(0)|
可見每次乘上這樣一個構造矩陣就能降次
可以推出：
(f(n),f(n-1),...,f(n-9))^(-1) = A^(n-9)*(f(9),f(8),...,f(0))^(-1)
?2.HDU 2294
好難。推好遞推公式后，難點就是怎么優化遞推過程，就是構造一個遞推矩陣，也要注意要經常取模
f[i][j]代表j種顏色組成i個珠子,有方程f[i][j]=f[i-1][j]*j+f[i-1][j-1]*(k-j-1),但由于N較大采用一般的方法DP必然按會超時的。
優化的方法是使用矩陣來做。將F[i-1]到F[i]的轉移用矩陣來描述，相當于一個k*k的線性變換矩陣。因此F[i]=A*F[i-1]，這里A是轉移矩陣，即F[i]=Ai-1*F[1]，所以F[1]+…+F[n]=A0*F[1]+…+An-1*F[1]=（E+A+A2+…+An-1)*F[1]。
0??k??0????0??0??0
0??1??k-1??0??0??0
0??0??2????0??0??0
.???.??.????.???.???.
.???.??.????.???.???.
0??0??0???0??k-1?1
0??0??0???0???0???k
for(i=1;i<=k;i++) //轉移矩陣A
? ? ? ? a[i-1][i]=k-i+1,?a[i][i]=i;
#include<cstdio> #include<string.h> #include <iostream> using namespace std; int cas, n, k; #define MAXN 40 #define LL long long const int mod = 1234567891; struct mat {int n, m;int data[MAXN][MAXN];mat() {n = m = 0;memset(data, 0, sizeof(data));} }; void debug(mat tmp) {for (int i = 0; i < tmp.n; i++) {for (int j = 0; j < tmp.m; j++)printf("%d ", tmp.data[i][j]);puts("");}puts(""); } mat mul(mat& a, mat& b) { //a乘b = cint i, j, k;mat c;c.n = a.n, c.m = b.m;for (i = 0; i < c.n; i++)for (j = 0; j < c.m; j++) {for (c.data[i][j] = k = 0; k < a.m; k++)c.data[i][j] = (c.data[i][j]+ ((LL) a.data[i][k] * b.data[k][j]) % mod) % mod;}return c; } mat pow(mat a, int x) {mat b;b.m = b.n = max(a.m, a.n);for (int i = 0; i < b.n; i++)b.data[i][i] = 1;while (x) {if (x & 1) {b = mul(a, b);}a = mul(a, a);x >>= 1;}return b; } mat add(mat a, mat b) { //矩陣乘法 mat c;c.m = a.m, c.n = a.n;for (int i = 0; i < a.n; i++)for (int j = 0; j < a.m; j++)c.data[i][j] = ((LL) a.data[i][j] + b.data[i][j]) % mod;return c; } //二分計算A^l+...A^r mat cal(mat A, int l, int r) {if (l == r)return pow(A, l);int m = (l + r + 1) / 2;mat t = cal(A, l, m - 1); //算前一半的和mat t2 = pow(A, m - l); //中介的矩陣t = add(t, mul(t2, t));if ((r - l + 1) & 1)t = add(t, pow(A, r)); //如果是奇數，則加上最后一個return t; }int main() { // freopen("data3.txt", "r", stdin);scanf("%d", &cas);while (cas--) {scanf("%d%d", &n, &k);mat A;A.m = A.n = k + 1;for (int i = 1; i <= k; i++) //轉移矩陣AA.data[i - 1][i] = k - i + 1, A.data[i][i] = i;printf("%d\n", cal(A, k, n).data[0][k]);}return 0; }?
3.高斯消元解線性方程
1.求有唯一解的線性方程組
#include<cstdio> #include<cmath> using namespace std; int n, m; double ans[55], d[55][55], eps = 1e-9; int gauss_cpivot(double a[55][55], int n, double x[]) {int i, j, k, p;for (j = 0; j < n; ++j) {for (i = j + 1, p = j; i < n; ++i)if (fabs(a[i][j]) > fabs(a[p][j]))p = i;if (fabs(a[p][j]) < eps)return 0;if (p != j)for (k = j; k <= n; ++k)swap(a[j][k], a[p][k]);for (i = j + 1; i < n; ++i)for (k = n; k >= j; --k)a[i][k] -= a[j][k] * a[i][j] / a[j][j];}for (j = n - 1; j >= 0; --j) {x[j] = a[j][n] / a[j][j];for (i = j - 1; i >= 0; --i)a[i][n] -= a[i][j] * x[j];}return 1; }?
2.求系數行列式的秩判斷是否有解
POJ 1830 用高斯消元解決開關問題:?每個開關只用一次，最后到達使若干個開關狀態改變
可以發現對于每個開關從起始狀態到終止狀態都是那些對它有影響的開關開啟關閉情況的異或，即用每個開關一個相關開關方程（相關為1，不相關為0）
那每個開關 系數以相關開關為1 建立一個異或的線性方程組 高斯消元看是否有解，有解每個解都是一種滿足條件的方法，則方法數就是自由未知量的所有狀態數
#include<cstdio> #include<cmath> #include<algorithm> using namespace std; #define MAXN 109 const double eps = 1e-8; inline bool ZERO(double a) {return fabs(a) < eps; } double d[55][55]; int R(double a[55][55], int r, int c) { //求系數行列式的秩int i, j, R, k, Max;for (R = j = 0; j < c - 1; j++) {for (i = R + 1, Max = R; i < r; i++)if (fabs(a[i][j]) > fabs(a[Max][j]))Max = i;if (Max != R)for (k = j; k < c; k++)swap(a[R][k], a[Max][k]);if (ZERO(a[R][j])) {continue;}for (i = R + 1; i < r; i++)if (!ZERO(a[i][j])) {double t = a[i][j] / a[R][j];for (k = j; k < c; k++)a[i][k] -= a[R][k] * t;}R++;}for (i = R; i < r; i++) { //無解情況if (!ZERO(a[i][c - 1]))return -1;}return R; } int start[MAXN]; int end[MAXN];int main() { // freopen("data3.txt", "r", stdin);int u, v;int T;int n;scanf("%d", &T);while (T--) {scanf("%d", &n);for (int i = 0; i < n; i++)scanf("%d", &start[i]);for (int i = 0; i < n; i++)scanf("%d", &end[i]);memset(d, 0, sizeof(d));while (scanf("%d%d", &u, &v) && u && v) {d[v - 1][u - 1] = 1;}for (int i = 0; i < n; i++)d[i][i] = 1;for (int i = 0; i < n; i++)d[i][n] = start[i] ^ end[i];int ans = R(d, n, n + 1);if (ans == -1)printf("Oh,it's impossible~!!\n");elseprintf("%d\n", 1 << (n - ans));}return 0; }?
http://www.matrix67.com/blog/archives/276
轉載于:https://www.cnblogs.com/updateofsimon/p/3409350.html
總結
以上是生活随笔為你收集整理的矩阵在计算机程序中的应用的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： python操作MongoDB
- 下一篇： java文件读写操作大全