📄 tmte.cpp
字号:
A(a+15,a+0) = 1; A(a+15,a+1) = -1;
A(a+16,a+14) = 1; A(a+16,a+0) = -1;
A(a+17,a+7) = 1; A(a+17,a+6) = -1;
A(a+18,a+13) = 1; A(a+18,a+5) = -1;
A(a+19,a+7) = 1; A(a+19,a+8) = -1;
A(a+20,a+17) = 1; A(a+20,a+3) = -1;
A(a+21,a+6) = 1; A(a+21,a+11) = -1;
Matrix U, V; DiagonalMatrix S;
SVD(A, S, U, V, true, true); CheckIsSorted(S);
DiagonalMatrix D(20); D = 1;
Matrix tmp = U.t() * U - D;
Clean(tmp,0.000000001); Print(tmp);
tmp = V.t() * V - D;
Clean(tmp,0.000000001); Print(tmp);
tmp = U * S * V.t() - A ;
Clean(tmp,0.000000001); Print(tmp);
}
{
Tracer et1("Stage 5");
Matrix A(10,10);
A.Row(1) << 1.00 << 0.07 << 0.05 << 0.00 << 0.06
<< 0.09 << 0.03 << 0.02 << 0.02 << -0.03;
A.Row(2) << 0.07 << 1.00 << 0.05 << 0.05 << -0.03
<< 0.07 << 0.00 << 0.07 << 0.00 << 0.02;
A.Row(3) << 0.05 << 0.05 << 1.00 << 0.05 << 0.02
<< 0.01 << -0.05 << 0.04 << 0.05 << -0.03;
A.Row(4) << 0.00 << 0.05 << 0.05 << 1.00 << -0.05
<< 0.04 << 0.01 << 0.02 << -0.05 << 0.00;
A.Row(5) << 0.06 << -0.03 << 0.02 << -0.05 << 1.00
<< -0.03 << 0.02 << -0.02 << 0.04 << 0.00;
A.Row(6) << 0.09 << 0.07 << 0.01 << 0.04 << -0.03
<< 1.00 << -0.06 << 0.08 << -0.02 << -0.10;
A.Row(7) << 0.03 << 0.00 << -0.05 << 0.01 << 0.02
<< -0.06 << 1.00 << 0.09 << 0.12 << -0.03;
A.Row(8) << 0.02 << 0.07 << 0.04 << 0.02 << -0.02
<< 0.08 << 0.09 << 1.00 << 0.00 << -0.02;
A.Row(9) << 0.02 << 0.00 << 0.05 << -0.05 << 0.04
<< -0.02 << 0.12 << 0.00 << 1.00 << 0.02;
A.Row(10) << -0.03 << 0.02 << -0.03 << 0.00 << 0.00
<< -0.10 << -0.03 << -0.02 << 0.02 << 1.00;
SymmetricMatrix AS; AS << A;
Matrix V; DiagonalMatrix D, D1;
ColumnVector Check(6);
EigenValues(AS,D,V); CheckIsSorted(D, true);
Check(1) = MaximumAbsoluteValue(A - V * D * V.t());
DiagonalMatrix I(10); I = 1;
Check(2) = MaximumAbsoluteValue(V * V.t() - I);
Check(3) = MaximumAbsoluteValue(V.t() * V - I);
EigenValues(AS, D1); CheckIsSorted(D1, true);
D -= D1;
Clean(D,0.000000001); Print(D);
Jacobi(AS,D,V);
Check(4) = MaximumAbsoluteValue(A - V * D * V.t());
Check(5) = MaximumAbsoluteValue(V * V.t() - I);
Check(6) = MaximumAbsoluteValue(V.t() * V - I);
SortAscending(D);
D -= D1;
Clean(D,0.000000001); Print(D);
Clean(Check,0.000000001); Print(Check);
// Check loading rows
SymmetricMatrix B(10);
B.Row(1) << 1.00;
B.Row(2) << 0.07 << 1.00;
B.Row(3) << 0.05 << 0.05 << 1.00;
B.Row(4) << 0.00 << 0.05 << 0.05 << 1.00;
B.Row(5) << 0.06 << -0.03 << 0.02 << -0.05 << 1.00;
B.Row(6) << 0.09 << 0.07 << 0.01 << 0.04 << -0.03
<< 1.00;
B.Row(7) << 0.03 << 0.00 << -0.05 << 0.01 << 0.02
<< -0.06 << 1.00;
B.Row(8) << 0.02 << 0.07 << 0.04 << 0.02 << -0.02
<< 0.08 << 0.09 << 1.00;
B.Row(9) << 0.02 << 0.00 << 0.05 << -0.05 << 0.04
<< -0.02 << 0.12 << 0.00 << 1.00;
B.Row(10) << -0.03 << 0.02 << -0.03 << 0.00 << 0.00
<< -0.10 << -0.03 << -0.02 << 0.02 << 1.00;
B -= AS; Print(B);
}
{
Tracer et1("Stage 6");
// badly scaled matrix
Matrix A(9,9);
A.Row(1) << 1.13324e+012 << 3.68788e+011 << 3.35163e+009
<< 3.50193e+011 << 1.25335e+011 << 1.02212e+009
<< 3.16602e+009 << 1.02418e+009 << 9.42959e+006;
A.Row(2) << 3.68788e+011 << 1.67128e+011 << 1.27449e+009
<< 1.25335e+011 << 6.05413e+010 << 4.34573e+008
<< 1.02418e+009 << 4.69192e+008 << 3.61098e+006;
A.Row(3) << 3.35163e+009 << 1.27449e+009 << 1.25571e+007
<< 1.02212e+009 << 4.34573e+008 << 3.69769e+006
<< 9.42959e+006 << 3.61098e+006 << 3.59450e+004;
A.Row(4) << 3.50193e+011 << 1.25335e+011 << 1.02212e+009
<< 1.43514e+011 << 5.42310e+010 << 4.15822e+008
<< 1.23068e+009 << 4.31545e+008 << 3.58714e+006;
A.Row(5) << 1.25335e+011 << 6.05413e+010 << 4.34573e+008
<< 5.42310e+010 << 2.76601e+010 << 1.89102e+008
<< 4.31545e+008 << 2.09778e+008 << 1.51083e+006;
A.Row(6) << 1.02212e+009 << 4.34573e+008 << 3.69769e+006
<< 4.15822e+008 << 1.89102e+008 << 1.47143e+006
<< 3.58714e+006 << 1.51083e+006 << 1.30165e+004;
A.Row(7) << 3.16602e+009 << 1.02418e+009 << 9.42959e+006
<< 1.23068e+009 << 4.31545e+008 << 3.58714e+006
<< 1.12335e+007 << 3.54778e+006 << 3.34311e+004;
A.Row(8) << 1.02418e+009 << 4.69192e+008 << 3.61098e+006
<< 4.31545e+008 << 2.09778e+008 << 1.51083e+006
<< 3.54778e+006 << 1.62552e+006 << 1.25885e+004;
A.Row(9) << 9.42959e+006 << 3.61098e+006 << 3.59450e+004
<< 3.58714e+006 << 1.51083e+006 << 1.30165e+004
<< 3.34311e+004 << 1.25885e+004 << 1.28000e+002;
SymmetricMatrix AS; AS << A;
Matrix V; DiagonalMatrix D, D1;
ColumnVector Check(6);
EigenValues(AS,D,V); CheckIsSorted(D, true);
Check(1) = MaximumAbsoluteValue(A - V * D * V.t()) / 100000;
DiagonalMatrix I(9); I = 1;
Check(2) = MaximumAbsoluteValue(V * V.t() - I);
Check(3) = MaximumAbsoluteValue(V.t() * V - I);
EigenValues(AS, D1);
D -= D1;
Clean(D,0.001); Print(D);
Jacobi(AS,D,V);
Check(4) = MaximumAbsoluteValue(A - V * D * V.t()) / 100000;
Check(5) = MaximumAbsoluteValue(V * V.t() - I);
Check(6) = MaximumAbsoluteValue(V.t() * V - I);
SortAscending(D);
D -= D1;
Clean(D,0.001); Print(D);
Clean(Check,0.0000001); Print(Check);
}
{
Tracer et1("Stage 7");
// matrix with all singular values close to 1
Matrix A(8,8);
A.Row(1)<<-0.4343<<-0.0445<<-0.4582<<-0.1612<<-0.3191<<-0.6784<<0.1068<<0;
A.Row(2)<<0.5791<<0.5517<<0.2575<<-0.1055<<-0.0437<<-0.5282<<0.0442<<0;
A.Row(3)<<0.5709<<-0.5179<<-0.3275<<0.2598<<-0.196<<-0.1451<<-0.4143<<0;
A.Row(4)<<0.2785<<-0.5258<<0.1251<<-0.4382<<0.0514<<-0.0446<<0.6586<<0;
A.Row(5)<<0.2654<<0.3736<<-0.7436<<-0.0122<<0.0376<<0.3465<<0.3397<<0;
A.Row(6)<<0.0173<<-0.0056<<-0.1903<<-0.7027<<0.4863<<-0.0199<<-0.4825<<0;
A.Row(7)<<0.0434<<0.0966<<0.1083<<-0.4576<<-0.7857<<0.3425<<-0.1818<<0;
A.Row(8)<<0.0<<0.0<<0.0<<0.0<<0.0<<0.0<<0.0<<-1.0;
Matrix U,V; DiagonalMatrix D;
SVD(A,D,U,V); CheckIsSorted(D);
Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
DiagonalMatrix I(8); I = 1; D -= I; Clean(D,0.0001); Print(D);
U *= U.t(); U -= I; Clean(U,0.000000001); Print(U);
V *= V.t(); V -= I; Clean(V,0.000000001); Print(V);
}
{
Tracer et1("Stage 8");
// check SortSV functions
Matrix A(15, 10);
int i, j;
for (i = 1; i <= 15; ++i) for (j = 1; j <= 10; ++j)
A(i, j) = i + j / 1000.0;
DiagonalMatrix D(10);
D << 0.2 << 0.5 << 0.1 << 0.7 << 0.8 << 0.3 << 0.4 << 0.7 << 0.9 << 0.6;
Matrix U = A; Matrix V = 10 - 2 * A;
Matrix Prod = U * D * V.t();
DiagonalMatrix D2 = D; SortDescending(D2);
DiagonalMatrix D1 = D; SortSV(D1, U, V); Matrix X = D1 - D2; Print(X);
X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
U = A; V = 10 - 2 * A;
D1 = D; SortSV(D1, U); X = D1 - D2; Print(X);
D1 = D; SortSV(D1, V); X = D1 - D2; Print(X);
X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
D2 = D; SortAscending(D2);
U = A; V = 10 - 2 * A;
D1 = D; SortSV(D1, U, V, true); X = D1 - D2; Print(X);
X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
U = A; V = 10 - 2 * A;
D1 = D; SortSV(D1, U, true); X = D1 - D2; Print(X);
D1 = D; SortSV(D1, V, true); X = D1 - D2; Print(X);
X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
}
{
Tracer et1("Stage 9");
// Tom William's example
Matrix A(10,10);
Matrix U;
Matrix V;
DiagonalMatrix Sigma;
Real myVals[] =
{
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
1, 1, 1, 1, 1, 1, 1, 1, 0, 0,
1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
};
A << myVals;
SVD(A, Sigma, U, V); CheckIsSorted(Sigma);
A -= U * Sigma * V.t();
Clean(A, 0.000000001); Print(A);
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -