问题描述
我试图计算任意大小的随机正交矩阵，但由于矩阵尺寸较小，因此面临着机器误差巨大的问题。我用Q^T * Q = I检查最终矩阵是否正交
其中Q是正交矩阵的计算值。例如，对于10 * 10矩阵的此操作将返回
1.000001421720586184 -0.000000728640227713 0.000001136830463799 -0.000000551609342727 -0.000001177027039965 0.000000334599582398 -0.000000858589995413 0.000000954985769303 0.000032744809653293 -0.000000265053286108
-0.000000728640227713 1.000000373167701527 -0.000000583888104495 0.000000285028920487 0.000000602479963850 -0.000000171504851561 0.000000439149502041 -0.000000489282575621 -0.000016836862737655 0.000000139458235281
0.000001136830463799 -0.000000583888104495 1.000000903071175539 -0.000000430043020645 -0.000000944743177025 0.000000267453533747 -0.000000690730534104 0.000000764348989692 0.000025922602192780 -0.000000194784469538
-0.000000551609342727 0.000000285028920487 -0.000000430043020645 1.000000193583126484 0.000000463290242271 -0.000000129639515781 0.000000340879278672 -0.000000371868992193 -0.000012221761904027 0.000000071060844115
-0.000001177027039965 0.000000602479963850 -0.000000944743177025 0.000000463290242271 1.000000972303993511 -0.000000277069934225 0.000000708304641621 -0.000000790186119274 -0.000027265457679301 0.000000229727452845
0.000000334599582398 -0.000000171504851561 0.000000267453533747 -0.000000129639515781 -0.000000277069934225 1.000000078745856667 -0.000000202136378957 0.000000224766235153 0.000007702164180343 -0.000000062098652538
-0.000000858589995413 0.000000439149502041 -0.000000690730534104 0.000000340879278672 0.000000708304641621 -0.000000202136378957 1.000000515565399593 -0.000000576213323968 -0.000019958173806416 0.000000172131276688
0.000000954985769303 -0.000000489282575621 0.000000764348989692 -0.000000371868992193 -0.000000790186119274 0.000000224766235153 -0.000000576213323968 1.000000641385961051 0.000022026878760393 -0.000000180133590665
0.000032744809653293 -0.000016836862737655 0.000025922602192780 -0.000012221761904027 -0.000027265457679301 0.000007702164180343 -0.000019958173806416 0.000022026878760393 1.000742765170839780 -0.000005353869978161
-0.000000265053286108 0.000000139458235281 -0.000000194784469538 0.000000071060844115 0.000000229727452845 -0.000000062098652538 0.000000172131276688 -0.000000180133590665 -0.000005353869978161 1.000000000000000000
所以我们可以看到矩阵是正交的，但非对角元素的误差很大 有什么解决办法吗？
我如何计算n * n个正交矩阵：
- 采用大小为1
Q = |1|的起始正方形正交矩阵 - 采用维度2
y = |rand(),rand()|的随机向量并将其规范为y = y/norm(y) - 用向量y构造一个Householder reflection并将其应用于右上角编号为1的矩阵Q，所以我有大小为2的正交矩阵Q。
- 在没有n * n矩阵的情况下重复执行，以增加维度的方式获取新的随机y。
代码：
#include <iostream>
#include <map>
#include <iostream>
#include <vector>
#include <iterator>
#include <cmath>
#include <utility>
using namespace std;
template<typename T>
T tolerance = T(1e-3);
template<typename T>
struct Triplet{
int i;
int j;
T b;
};
template<typename T>
T Tabs(T num){
if(num<T(0)) return -num;
else return num;
}
template<typename T>
class DOK{
private:
/*
* Dictionary of Keys,pair<int,int> is coordinates of non-zero elements,* next int is value
*/
int size_n;
int size_m;
map<pair<int,int>,T> dict;
// int count;
public:
DOK(vector<Triplet<T>> &matrix,int n,int m){
this->resize(n,m);
this->fill(matrix);
}
DOK(int n,m);
}
~DOK() = default;
void fill(vector<Triplet<T>> &matrix){
//this->count=matrix.size();
//cout<<"Input your coordinates with value in format \"i j val\" "<<endl;
for(int k = 0; k < matrix.size(); k++){
this->insert(matrix[k]);
}
}
void insert(const Triplet<T> &Element){
if(Element.i >= this->size_n){
this->size_n = Element.i+1;
}
if(Element.j >= this->size_m){
this->size_m = Element.j+1;
}
pair<int,int> coordinates = {Element.i,Element.j};
this->dict.insert(pair(coordinates,Element.b));
}
void resize(int n,int m){
this->size_n=n;
this->size_m=m;
}
void print() const{
cout<<endl;
for(int i = 0; i < this->size_n; i++){
for(int j = 0; j < this->size_m; j++){
if(this->dict.find({i,j})!= this->dict.cend()) cout<< this->dict.find(pair(i,j))->second<<" "; else cout<<0<<" ";
}
cout<<endl;
}
}
void clearZeros(){
for(auto i = this->dict.begin(); i!=this->dict.end();){
if(Tabs(i->second) <= tolerance<T>){
i = this->dict.erase(i);
} else{
i++;
}
}
}
[[nodiscard]] pair<int,int> getSize() const{
return {size_n,size_m};
}
DOK<T> transpose(){
DOK<T> A = DOK<T>(this->size_m,this->size_n);
for(auto &i: this->dict){
A.insert({i.first.second,i.first.first,i.second});
}
return A;
}
DOK<T>& operator-=(const DOK<T>& matrix){
try{
if(this->size_n != matrix.size_n || this->size_m != matrix.size_m) throw 1;
for(auto j: matrix.dict){
if(this->dict.find(j.first)!=this->dict.cend()) {
this->dict[j.first] -= j.second;
}else{
this->dict.insert({j.first,-j.second});
//M.count++;
}
}
this->clearZeros();
return *this;
}
catch (int a) {
cout<<"Sizes of Matrices are different."<<endl;
}
}
DOK<T> operator-(const DOK<T> &matrix) const{
DOK<T> t = *this;
return move(t-=matrix);
}
DOK<T>& operator*=(const DOK<T> &matrix){
try {
if(this->size_m != matrix.size_n) throw 1;
DOK<T> M = DOK(this->size_n,matrix.size_m);
for (int i = 0; i < this->size_n; i++) {
for (int j = 0; j < matrix.size_m; j++) {
T a=0;
for(int k = 0; k<this->size_m; k++){
if(this->dict.find({i,k}) != this->dict.cend() && matrix.dict.find({k,j})!=matrix.dict.cend()){
a+=this->dict.find({i,k})->second*matrix.dict.find({k,j})->second;
//cout<<a<<endl;
}
}
Triplet<T> m = {i,j,a};
M.insert(m);
}
}
this->clearZeros();
*this=M;
return *this;
}
catch (int a) {
cout<<"Wrong sizes of matrices to multiplication"<<endl;
}
}
DOK<T> operator*(const DOK<T>& matrix) const{
DOK<T> t = *this;
return t*=matrix;
}
DOK<T>& operator*=(T& k){
for(auto i: this->dict){
this->dict[i.first]*=k;
}
this->clearZeros();
return *this;
}
DOK<T> operator*(T& k) const{
DOK<T> t = *this;
return move(t*=k);
}
DOK<T>& operator*=(const T& k){
for(auto i: this->dict){
this->dict[i.first]*=k;
}
this->clearZeros();
return *this;
}
};
template<typename T>
vector<T> operator*(const DOK<T> &matrix,const vector<T> &x){
vector<T> result;
for(int i = 0; i < x.size(); i++){
T temp = 0;
for(int j = 0; j < x.size(); j++){
temp+=matrix(i,j)*x[j];
}
result.push_back(temp);
}
return move(result);
}
template<typename T>
T operator*(const vector<T> &x,const vector<T> &b) {
T result = 0;
for(int i = 0; i < x.size(); i++){
result+=x[i]*b[i];
}
}
template<typename T>
vector<T> operator*(const vector<T> &x,const DOK<T> &matrix) {
vector<T> result;
for(int i = 0; i < x.size(); i++){
T temp = 0;
for(int j = 0; j < x.size(); j++){
temp+=matrix(j,i)*x[j];
}
result.push_back(temp);
}
return move(result);
}
template<typename T>
DOK<T> operator*(T& k,const DOK<T> &matrix){
return matrix*k;
}
template<typename T>
DOK<T> operator*(const T& k,const DOK<T> &matrix){
return matrix*k;
}
template<typename T>
vector<T>& operator*=(const DOK<T> &matrix,const vector<T> &x){
return matrix*x;
}
template<typename T>
vector<T>& operator*=(const vector<T> &x,const DOK<T> &matrix){
return x*matrix;
}
template<typename T>
vector<T> operator*(const vector<T> &x,T k){
vector<T> result = x;
for(int i = 0; i<x.size(); i++){
result[i]*=k;
}
return result;
}
template<typename T>
vector<T> operator*(T k,const vector<T> &x) {
return x*k;
}
template<typename T>
ostream& operator<<(ostream &os,const DOK<T> &matrix) {
matrix.DOK<T>::print();
return os;
}
template<typename T>
T norm(const vector<T> x){
T result = 0;
for(int i = 0; i < x.size(); i++){
result+=pow(x[i],2);
}
return sqrt(result);
}
template<typename T>
DOK<T> Ortogonal(int n){
srand(time(NULL));
vector<Triplet<T>> in = {{0,T(1)}};
DOK<T> Q = DOK<T>(in,1,1);
DOK<T> E = Q;
vector<T> y;
for(int i = 1; i<n; i++){
y.clear();
for(int m = 0; m<i+1; m++){
y.emplace_back(rand());
}
y = (1/norm(y))*y;
DOK<T> Y = DOK<T>(i+1,i+1);
for(int j = 0; j<i+1; j++){
for(int k = 0; k<i+1; k++){
Y.insert({j,k,y[j]*y[k]});
}
}
Q.insert({i,i,T(1)});
cout<<Q;
Y*=T(2);
E.insert({i,T(1)});
Q = (E - Y)*Q;
}
return Q;
}
main.cpp：
#include <iostream>
#include "DOK.h"
using namespace std;
int main() {
DOK<long double> O = Ortogonal<long double>(10);
cout<<O.transpose()*O;
return 0;
}
DOK是具有所有重载运算符的稀疏矩阵的模板类。
解决方法
我试图重现该问题，并且似乎在基于Householder的过程中，程序每次执行时的数字错误量都会发生变化。
这意味着由于某种原因，该错误对所使用的随机数序列很敏感。该顺序取决于通过srand(time(NULL))函数调用的程序的确切启动时间。
如果“住户程序”不是必要条件，我建议改用Gram-Schmidt based procedure。
该过程如下：
1) create a full n*n random square matrix
2) normalize its column vectors
3) for each column vector except the leftmost one:
substract from that vector its projections on column vectors located to its left
renormalize the column vector
不幸的是，我对DOK稀疏矩阵框架不够熟悉，无法在其中解决这样的问题。但这当然可以由更熟悉它的程序员来完成。
改为使用Eigen开源线性代数库，代码如下所示：
#include <vector>
#include <random>
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
MatrixXd makeRandomOrthogonalMatrix(int dim,int seed)
{
std::normal_distribution<double> dist(0.0,1.0);
std::mt19937 gen(seed);
std::function<double(void)> rng = [&dist,&gen] () { return dist(gen); };
MatrixXd rawMat = MatrixXd::NullaryExpr(dim,dim,rng);
// massage rawMat into an orthogonal matrix - Gram-Schmidt process:
for (int j=0; j < dim; j++) {
VectorXd colj = rawMat.col(j);
colj.normalize();
rawMat.col(j) = colj;
}
for (int j=1; j < dim; j++) {
VectorXd colj = rawMat.col(j);
for (int k = 0; k < j; k++) {
VectorXd colk = rawMat.col(k);
colj -= (colk.dot(colj)) * colk;
}
colj.normalize();
rawMat.col(j) = colj;
}
return rawMat;
}
关于伪随机数产生的注释：
DOK程序中的函数rand()在零和大整数之间提供均匀的分布。在这里，由于以下原因，我切换到了高斯（也称为“正态”）分布：
- 这是相关数学文献所使用的
- 这确保随机列向量没有特权方向
比率对称性仅来自表达式exp（-x 2 ）* exp（-y 2 ）* exp（-z 2 ）等于exp（-（x 2 + y 2 + z 2 ）））。
使用的测试程序：
int main()
{
int ranSeed = 4241;
const int dim = 10;
MatrixXd orthoMat = makeRandomOrthogonalMatrix(dim,ranSeed);
MatrixXd trOrthoMat = orthoMat.transpose();
MatrixXd productMat2 = trOrthoMat * orthoMat;
std::cout << "productMat2 = \n" << productMat2 << '\n' << std::endl;
return EXIT_SUCCESS;
}
即使为随机种子尝试了多个值，数值误差的数量也看起来令人满意：
productMat2 =
1 1.38778e-17 0 1.38778e-17 -9.71445e-17 1.79544e-16 1.11022e-16 1.80411e-16 4.16334e-17 -2.98372e-16
1.38778e-17 1 0 -6.93889e-17 -1.38778e-17 5.89806e-17 5.55112e-17 -4.16334e-17 -2.77556e-17 -7.63278e-17
0 0 1 0 -1.38778e-17 7.37257e-17 0 9.71445e-17 1.80411e-16 -1.38778e-17
1.38778e-17 -6.93889e-17 0 1 1.38778e-17 -1.64799e-17 -5.55112e-17 0 -5.55112e-17 2.77556e-17
-9.71445e-17 -1.38778e-17 -1.38778e-17 1.38778e-17 1 -4.25007e-17 5.55112e-17 -1.38778e-17 -6.93889e-17 1.78677e-16
1.79544e-16 5.89806e-17 7.37257e-17 -1.64799e-17 -4.25007e-17 1 7.80626e-17 -5.37764e-17 2.60209e-17 -9.88792e-17
1.11022e-16 5.55112e-17 0 -5.55112e-17 5.55112e-17 7.80626e-17 1 8.32667e-17 8.32667e-17 -7.63278e-17
1.80411e-16 -4.16334e-17 9.71445e-17 0 -1.38778e-17 -5.37764e-17 8.32667e-17 1 4.16334e-17 1.04083e-16
4.16334e-17 -2.77556e-17 1.80411e-16 -5.55112e-17 -6.93889e-17 2.60209e-17 8.32667e-17 4.16334e-17 1 -7.28584e-17
-2.98372e-16 -7.63278e-17 -1.38778e-17 2.77556e-17 1.78677e-16 -9.88792e-17 -7.63278e-17 1.04083e-16 -7.28584e-17 1
在非特征环境中使用该功能：
要允许基于DOK的程序使用该随机unit矩阵，我们可以使用常规的std::vector对象，以使DOK类型系统与本征类型系统隔离。使用以下管道代码：
// for export to non-Eigen programs:
std::vector<double> makeRandomOrthogonalMatrixAsVector(int dim,int seed)
{
MatrixXd randMat = makeRandomOrthogonalMatrix(dim,seed);
std::vector<double> matVec;
for (int i=0; i < dim; i++)
for (int j=0; j < dim; j++)
matVec.push_back(randMat(i,j));
return matVec;
}
最后一步，可以在事物的DOK方面重建随机unit矩阵，而不会过度暴露于特征类型。此代码可用于从中间向量重建矩阵：
extern vector<double> makeRandomOrthogonalMatrixAsVector(int,int);
int main(int argc,const char* argv[])
{
int dim = 10;
int randomSeed = 4241;
vector<double> vec = makeRandomOrthogonalMatrixAsVector(dim,randomSeed);
// create matrix in DoK format:
vector<Triplet<double>> trv;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++) {
double x = vec[i*dim + j];
trv.push_back(Triplet<double>{i,j,x});
}
DOK<double> orthoMat(trv,dim);
DOK<double> prodMat = orthoMat.transpose() * orthoMat;
cout << "\nAlmost Unit Matrix = \n" << prodMat << std::endl;
return 0;
}
DOK程序输出：
Almost Unit Matrix =
1 0 4.16334e-17 1.11022e-16 -8.32667e-17 -2.77556e-17 0 -2.77556e-17 -1.38778e-17 8.32667e-17
0 1 -1.5786e-16 -1.38778e-16 -5.55112e-17 2.42861e-17 1.11022e-16 2.22045e-16 -5.96745e-16 -2.77556e-17
4.16334e-17 -1.5786e-16 1 1.50487e-16 6.93889e-18 5.55112e-17 -5.72459e-17 -8.32667e-17 5.10009e-16 -1.83881e-16
1.11022e-16 -1.38778e-16 1.50487e-16 1 -6.93889e-18 -1.21431e-17 1.73472e-17 0 1.45717e-16 2.94903e-17
-8.32667e-17 -5.55112e-17 6.93889e-18 -6.93889e-18 1 -9.36751e-17 -2.77556e-16 -2.498e-16 5.68989e-16 -2.91434e-16
-2.77556e-17 2.42861e-17 5.55112e-17 -1.21431e-17 -9.36751e-17 1 1.04083e-16 -3.46945e-17 1.04083e-17 -3.40006e-16
0 1.11022e-16 -5.72459e-17 1.73472e-17 -2.77556e-16 1.04083e-16 1 -1.38778e-16 2.77556e-16 -1.04083e-16
-2.77556e-17 2.22045e-16 -8.32667e-17 0 -2.498e-16 -3.46945e-17 -1.38778e-16 1 0 9.71445e-17
-1.38778e-17 -5.96745e-16 5.10009e-16 1.45717e-16 5.68989e-16 1.04083e-17 2.77556e-16 0 1 2.08167e-17
8.32667e-17 -2.77556e-17 -1.83881e-16 2.94903e-17 -2.91434e-16 -3.40006e-16 -1.04083e-16 9.71445e-17 2.08167e-17 1
因此，正如预期的那样，数值误差仍然很小。
请注意，与您发布的DOK程序不同，除非在源代码中手动更改了随机种子的值，否则使用的随机数序列始终保持不变。一个常见的改进是使用int seed = std::stoi(argv[1]);之类的代码将随机种子作为命令行参数传递。