0

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!

Michael E2
  • 235,386
  • 17
  • 334
  • 747
John hall
  • 31
  • 4

0 Answers0