1
$\begingroup$

I am trying to write Mathematica code to find $R $ (upper triangular matrix) in $ A=QR $ $A \in \Bbb{R}^{n\times n}$decomposition using reflection method. I wrote for $3 \times 3$ matrix. Any tips for write it for $n$ dimension.

$\endgroup$
4
  • $\begingroup$ How did you implement your Householder reflections? $\endgroup$ Dec 20, 2016 at 10:18
  • $\begingroup$ Finding the rotation angle and Q[b_] := Module[{[Theta]}, [Theta] = Theta[b]; R[[Theta]].( { {1, 0}, {0, -1} } ).R[-[Theta]]] $\endgroup$
    – Leonard
    Dec 20, 2016 at 10:21
  • $\begingroup$ Why not try to write a method for the Householder reflection for $n$-dimensions before anything else? $\endgroup$ Dec 20, 2016 at 10:23
  • $\begingroup$ Have you see the Method option of Orthogonalize and the Matrix Decompositions guide? $\endgroup$
    – Edmund
    Dec 20, 2016 at 13:55

1 Answer 1

5
$\begingroup$

Here's a pretty basic (and inefficient) routine for getting the QR decomposition of a matrix via Householder reflection:

qrd[mat_?MatrixQ] := Module[{r = mat, h, m, n, q, v, v2},
    {m, n} = Dimensions[r]; q = IdentityMatrix[m];
    Do[v = PadLeft[r[[k ;;, k]], m];
       v2 = v - SparseArray[{k -> Norm[v]}, m];
       h = If[! TrueQ[Norm[v2, ∞] == 0], ReflectionMatrix[v2], IdentityMatrix[m]];
       q = q.h; r = h.r,
       {k, n}];
    {q, r}]

which returns the $\mathbf Q$ and $\mathbf R$ factors in a list.

It could be made more efficient by exploiting the structure of the Householder update (see e.g. Golub and Van Loan or Stewart for how to do that), but I'll leave that for the OP or someone else to do.

$\endgroup$
2
  • $\begingroup$ Thank you and I am trying to do so. $\endgroup$
    – Leonard
    Dec 20, 2016 at 12:20
  • 2
    $\begingroup$ Of course, since QRDecomposition[] is actually built-in, the code here is more of a toy than anything else. $\endgroup$ Dec 20, 2016 at 14:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.