So I have this code for QR Householder decomposition:
QRDecHouseholder[matrix_, size_] :=
Module[{i, j, Avector1, vectoridentity, Sumcolumnvector, vectork, hk,
H, Q, R},
R = matrix;
H = {};
For[i = 1, i <= size[[2]] - 1, i++,
Avector1 = R[[i ;; size[[2]], i]];
vectoridentity = IdentityMatrix[Length[Avector1]][[1]];
Sumcolumnvector = Total[(Avector1)^2];
vectork =
Sign[Avector1[[1]]] * (Sqrt[Sumcolumnvector]) * vectoridentity +
Avector1 ;
hk = (IdentityMatrix[Length[Avector1]]) -
2*((Outer[Times, vectork, vectork])/(Dot[vectork, vectork]));
If[i < 2, ,
hk = ArrayFlatten[{{IdentityMatrix[i - 1],
ConstantArray[0, {i - 1, size[[2]] - i + 1}]}, {ConstantArray[
0, {size[[2]] - i + 1, i - 1}], hk}}]];
R = Dot[hk, R];
AppendTo[H, hk];
];
Q = ParallelCombine[Dot[Sequence @@ ##] &, H, Dot];
Return[List[Q, Chop[R]]]
]
And I think it's working perfectly. It recieves a matrix and the Dimensions[matrix] which is the size.
Now I have this code to find Q transposed, An (upper triangular matrix) and Q:
NewDecompMTM[matrix_, size_, itermax_, abserr_] :=
Module[{i, Aiter, Aiterb4, Q, Qiter, R, valerr},
Aiter = matrix;
valerr = 1000;
Q = IdentityMatrix[size[[2]]];
For[i = 1, valerr > abserr && i < itermax, i++,
Aiterb4 = Aiter;
{Qiter, R} = QRDecHouseholder[Aiter, size];
Aiter = Dot[R, Qiter];
Q = Dot[Transpose[Qiter], Q];
valerr = Max[Abs[Diagonal[Aiter] - Diagonal[Aiterb4]]];
];
Return[List[Transpose[Q], Chop[Aiter], Q]]
]
But when I try with a simple matrix like
A = {{1, 2, 3}, {4, 5, 6}, {3, 4, 6}};
Dimensions[A]
NewDecompMTM[A, Dimensions[A], 30, 0.0005]
Householder works but NewDecMTM doesn't and gives the error:
N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -(72/(25+(1+Sqrt[26])^2))-(8 (1+Sqrt[26]))/(25+(1+Sqrt[26])^2)+4 (1-32/(25+Plus[<<2>>]^2)).
N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -(72/(25+(1+Sqrt[26])^2))-(8 (1+Sqrt[26]))/(25+(1+Sqrt[26])^2)+4 (1-32/(25+Plus[<<2>>]^2)).
N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -(96/(25+(1+Sqrt[26])^2))-(6 (1+Sqrt[26]))/(25+(1+Sqrt[26])^2)+3 (1-18/(25+Plus[<<2>>]^2)).
General::stop: Further output of N::meprec will be suppressed during this calculation.
I need your help. Thanks!
SingularValueDecomposition? – flinty Dec 06 '20 at 14:54