游戏网站建设公司,建设银行网站登陆二星是什么意思,可以免费学编程的网站,wordpress编辑富文任务
生成一个4 3的随机矩阵A#xff0c;应用Jacobi方法求解矩阵 A A T \mathbf{AA^T} AAT的特征值#xff0c;计算矩阵A的SVD分解。要求A的每个元素均为[0; 1]区间内的随机数。
算法
Jacobi方法求矩阵 A A T \mathbf{AA^T} AAT的特征值;将 A A T \mathbf{AA^T} AAT的特征…任务
生成一个4 × 3的随机矩阵A应用Jacobi方法求解矩阵 A A T \mathbf{AA^T} AAT的特征值计算矩阵A的SVD分解。要求A的每个元素均为[0; 1]区间内的随机数。
算法
Jacobi方法求矩阵 A A T \mathbf{AA^T} AAT的特征值;将 A A T \mathbf{AA^T} AAT的特征值从大到小排序后用高斯消元法求解每个特征值对应的特征向量并归一化然后转置为列向量组得到 U 4 × 4 \mathbf{U}_{4\times4} U4×4;矩阵 Σ 4 × 3 \mathbf{\Sigma}_{4\times3} Σ4×3的主对角元为 A A T \mathbf{AA^T} AAT的非零特征值的算术平方根本题中仅考虑了实数情况其余位置全为0;矩阵 V T 3 × 3 \mathbf{V^T}_{3\times3} VT3×3每行即 U T \mathbf{U^T} UT A \mathbf{A} A对应行除以 Σ \mathbf{\Sigma} Σ对应的主对角元。
注理论上 A A T \mathbf{AA^T} AAT和 A T A \mathbf{A^TA} ATA的非零特征值相同但是程序计算出的结果会有一些差异因此 Σ \mathbf{\Sigma} Σ的主对角元不能取 A T A \mathbf{A^TA} ATA的特征值的算术平方根 U T \mathbf{U^T} UT的每行元素也不能直接取 A T A \mathbf{A^TA} ATA的特征向量而要按第4步所述计算详情请参考这篇论文。
代码
主函数
int main()
{vectorvectorlong double A generateRandomMatrix(4, 3);cout A: endl;for(int i 0; i 4; i){for(int j 0; j 3; j)cout A[i][j] ;cout endl;}cout endl;vectorvectorlong double AT calAT(A);vectorvectorlong double AAT multiplyMatrices(A, AT);int n1 AAT.size();int n2 A[0].size();vectorlong double x calEigenValue(AAT);sort(x.begin(), x.end());reverse(x.begin(), x.end());vectorvectorlong double U calAT(Normalization(calEigenVector(AAT, x)));cout U: endl;for(int i 0; i n1; i){for(int j 0; j n1; j)cout U[i][j] ;cout endl;}cout endl;vectorvectorlong double Sigma1 calSigma(x, n1, n2);cout Sigma: endl;for(int i 0; i n1; i){for(int j 0; j n2; j)cout Sigma1[i][j] ;cout endl;}cout endl;vectorvectorlong double VT calculateVT(U, A, Sigma1, n2);cout VT: endl;for(int i 0; i n2; i){for(int j 0; j n2; j)cout VT[i][j] ;cout endl;}return 0;
}生成一个具有指定行数和列数的随机矩阵
vectorvectorlong double generateRandomMatrix(int rows, int cols) {random_device rd;//随机数种子mt19937 gen(rd());//使用 rd 生成的随机数来初始化一个mt19937类型的随机数生成器 genuniform_real_distribution dis(0.0, 1.0);//创建一个均匀实数分布 disvectorvectorlong double matrix(rows, vectorlong double(cols));for (int i 0; i rows; i)for (int j 0; j cols; j)matrix[i][j] dis(gen);//生成随机数return matrix;
}找到矩阵 A 中非对角元中绝对值最大的元素并返回其位置
pairint, int chooseMax(vectorvectorlong double A)
{long double max 0;pairint, int maxPos;int n A.size();for(int i 0; i n; i)for(int j 0; j n; j)if(i ! j fabsl(A[i][j]) max){max fabsl(A[i][j]);maxPos make_pair(i, j);}return maxPos;
}计算矩阵 A 的转置
vectorvectorlong double calAT(vectorvectorlong double A)
{int n1 A.size();int n2 A[0].size();vectorvectorlong double AT(n2, vectorlong double(n1));for(int i 0; i n1; i)for(int j 0; j n2; j)AT[j][i] A[i][j];return AT;
}计算两个矩阵的乘积
vectorvectorlong double multiplyMatrices(const vectorvectorlong double A, const vectorvectorlong double B) {int n1 A.size();int n2 B[0].size();int n3 A[0].size();vectorvectorlong double result(n1, vectorlong double(n2, 0.0));for(int i 0; i n1; i) {for(int j 0; j n2; j) {for(int k 0; k n3; k) {result[i][j] A[i][k] * B[k][j];}}}return result;
}计算矩阵 Q^T * A * Q
vectorvectorlong double calQTAQ(vectorvectorlong double A)
{int n A.size();pairint, int maxPos chooseMax(A);int row maxPos.first;int col maxPos.second;long double s (A[col][col] - A[row][row]) / (2 * A[row][col]);long double t 0;if(s 0)t 1;else if(abs(-s sqrt(1 s * s)) abs(-s - sqrt(1 s * s)))t -s sqrt(1 s * s);elset -s - sqrt(1 s * s);long double c 1 / sqrt(1 t * t);long double d t * c;vectorvectorlong double QTAQ(n, vectorlong double(n));for(int i 0; i n; i)for(int j 0; j n; j)QTAQ[i][j] calculateElement(A, i, j, row, col, t, c, d);return QTAQ;
}// 计算矩阵Q^T * A * Q的每个元素使用给定的参数 p, q, t, c, d
long double calculateElement(const vectorvectorlong double A, int i, int j, long double p, long double q, long double t, long double c, long double d) {if (i p j p)return A[p][p] - t * A[p][q];else if (i q j q)return A[q][q] t * A[p][q];else if ((i p j q) || (i q j p))return 0;else if (i ! q i ! p (j p || j q))return (j p ? c : d) * A[p][i] - (j p ? d : (-c)) * A[q][i];else if ((i p || i q) j ! q j ! p)return (i p ? c : d) * A[p][j] - (i p ? d : (-c)) * A[q][j];elsereturn A[i][j];
}判断Jacobi方法是否满足结束条件
int judgeEnd(vectorvectorlong double A)
{int i, j;int n A.size();for(i 0; i n; i)for(j 0; j n; j)if(i ! j fabsl(A[i][j]) 1e-6)return 0;if(i n j n) return 1;
}Jacobi方法计算矩阵 A 的特征值
vectorlong double calEigenValue(vectorvectorlong double A)
{int n A.size();vectorlong double eigenValue(n);vectorvectorlong double QTAQ calQTAQ(A);int i, j;while(!judgeEnd(QTAQ))QTAQ calQTAQ(QTAQ);for(i 0; i n; i)eigenValue[i] QTAQ[i][i];return eigenValue;
}使用给定的特征值计算矩阵 A 的特征向量
vectorvectorlong double calEigenVector(vectorvectorlong double A, vectorlong double eigenValue)
{int n A.size();int num 0;vectorvectorlong double x(n, vectorlong double(n));vectorvectorlong double tempMartix(n, vectorlong double(n));vectorvectorlong double eigenVector(n, vectorlong double(n));for(int k 0; k n; k){for(int i 0; i n; i)for(int j 0; j n; j)i j ? tempMartix[i][j] A[i][j] - eigenValue[k] : tempMartix[i][j] A[i][j];vectorvectorlong double B Column_Elimination(tempMartix);int cnt 0;//记录消元后全为0的行数for(int i 0; i n; i){for(int j 0; j n; j){if(fabsl(B[i][j]) 1e-7)break;else if(j n - 1)cnt;}}vectorvectorlong double result solve(B, cnt);for(int i 0; i cnt; i)copy(result[i].begin(), result[i].end(), x[num i].begin());num cnt;}return x;
}求解有无穷多解的线性方程组
// 对矩阵A进行列主元化成上三角
vectorvectorlong double Column_Elimination(vectorvectorlong double A)
{int n A.size();vectorvectorlong double Temp(n, vectorlong double(n));for(int i 0; i n; i)for(int j 0; j n; j)Temp[i][j] A[i][j];for(int col 0; col n; col){long double maxnum abs(Temp[col][col]);int maxrow col;for(int row col 1; row n; row){if(abs(Temp[row][col]) maxnum){maxnum abs(Temp[row][col]);maxrow row;}}swap(Temp[col], Temp[maxrow]);for(int row col 1; row n; row){long double res Temp[row][col] / Temp[col][col];for(int loc col; loc n; loc)Temp[row][loc] - Temp[col][loc] * res; }}return Temp;
}// 求解系数矩阵为上三角矩阵A且A的行列式不为0的线性方程组
vectorlong double SolveUpperTriangle(vectorvectorlong double A, vectorlong double b)
{int n A.size();vectorlong double x(n);vectorvectorlong double Temp(n, vectorlong double(n1));for(int i 0; i n; i)for(int j 0; j n; j)Temp[i][j] A[i][j];for(int i 0; i n; i)Temp[i][n] b[i];for(int row n-1; row 0; row--){for(int col row 1; col n; col){Temp[row][n] - Temp[col][n] * Temp[row][col] / Temp[col][col];Temp[row][col] 0;}Temp[row][n] / Temp[row][row];Temp[row][row] 1;}for(int i 0; i n; i)x[i] Temp[i][n];return x;
}// 解系数矩阵为上三角矩阵 A 的线性方程组且A全为0的行数为 cnt
vectorvectorlong double solve(vectorvectorlong double A, int cnt)
{int n A.size();vectorvectorlong double x(cnt, vectorlong double(n));vectorvectorlong double Temp(n-cnt, vectorlong double(n-cnt));vectorlong double Tempb(n-cnt);for(int i 0; i cnt; i){for(int j n - 1; j n - cnt; j--){if(j n - i)x[i][j] 0;elsex[i][j] 1;}}for(int i 0; i n - cnt; i)for(int j 0; j n - cnt; j)Temp[i][j] A[i][j];for(int i 0; i cnt; i){for(int j n - cnt - 1; j 0; j--){Tempb[j] 0;for(int k 0; k cnt; k)Tempb[j] - A[j][n- cnt k] * x[i][n- cnt k];}vectorlong double res SolveUpperTriangle(Temp, Tempb);for(int j 0; j n - cnt; j)x[i][j] res[j];}return x;
}使用给定的特征值 x 和矩阵的行数 n1 和列数 n2计算 Sigma 矩阵
vectorvectorlong double calSigma(vectorlong double x, int n1, int n2)
{vectorvectorlong double Sigma(n1, vectorlong double(n2));for(int i 0; i min(n1, n2); i)Sigma[i][i] sqrt(x[i]);return Sigma;
}计算向量 x 的欧几里得范数
long double EuclideanNorm(vectorlong double x)
{long double norm 0;for(int i 0; i x.size(); i)norm x[i] * x[i];return sqrt(norm);
}对矩阵 A 进行归一化
vectorvectorlong double Normalization(vectorvectorlong double A)
{int rows A.size();for(int i 0; i rows; i){long double norm EuclideanNorm(A[i]);int cols A[i].size();for(int j 0; j cols; j)A[i][j] / norm;}return A;
}计算 VT 矩阵
vectorvectorlong double calculateVT(const vectorvectorlong double U, const vectorvectorlong double A, const vectorvectorlong double Sigma1, int n2) {vectorvectorlong double UTA multiplyMatrices(calAT(U), A);vectorvectorlong double VT(n2, vectorlong double(n2));for(int i 0; i n2; i)for(int j 0; j n2; j)VT[i][j] UTA[i][j]/Sigma1[i][i];return VT;
}
