Problem Statement

Any rectangular

m × n

{\displaystyle m\times n}

matrix

A

{\displaystyle {\mathit {A}}}

can be decomposed to a product of an orthogonal matrix

Q

{\displaystyle {\mathit {Q}}}

and an upper (right) triangular matrix

R

{\displaystyle {\mathit {R}}}

, as described in QR decomposition. Task Demonstrate the QR decomposition on the example matrix from the Wikipedia article: and the usage for linear least squares problems on the example from Polynomial regression. The method of Householder reflections should be used: Method Multiplying a given vector

a

{\displaystyle {\mathit {a}}}

, for example the first column of matrix

A

{\displaystyle {\mathit {A}}}

, with the Householder matrix

H

{\displaystyle {\mathit {H}}}

, which is given as reflects

a

{\displaystyle {\mathit {a}}}

about a plane given by its normal vector

u

{\displaystyle {\mathit {u}}}

. When the normal vector of the plane

u

{\displaystyle {\mathit {u}}}

is given as then the transformation reflects

a

{\displaystyle {\mathit {a}}}

onto the first standard basis vector which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of

a

1

{\displaystyle a_{1}}

: and normalize with respect to the first element: The equation for

H

{\displaystyle H}

thus becomes: or, in another form with Applying

H

{\displaystyle {\mathit {H}}}

on

a

{\displaystyle {\mathit {a}}}

then gives and applying

H

{\displaystyle {\mathit {H}}}

on the matrix

A

{\displaystyle {\mathit {A}}}

zeroes all subdiagonal elements of the first column: In the second step, the second column of

A

{\displaystyle {\mathit {A}}}

, we want to zero all elements but the first two, which means that we have to calculate

H

{\displaystyle {\mathit {H}}}

with the first column of the submatrix (denoted *), not on the whole second column of

A

{\displaystyle {\mathit {A}}}

. To get

H

2

{\displaystyle H_{2}}

, we then embed the new

H

{\displaystyle {\mathit {H}}}

into an

m × n

{\displaystyle m\times n}

identity: This is how we can, column by column, remove all subdiagonal elements of

A

{\displaystyle {\mathit {A}}}

and thus transform it into

R

{\displaystyle {\mathit {R}}}

. The product of all the Householder matrices

H

{\displaystyle {\mathit {H}}}

, for every column, in reverse order, will then yield the orthogonal matrix

Q

{\displaystyle {\mathit {Q}}}

. The QR decomposition should then be used to solve linear least squares (Multiple regression) problems

A

x

b

{\displaystyle {\mathit {A}}x=b}

by solving When

R

{\displaystyle {\mathit {R}}}

is not square, i.e.

m

n

{\displaystyle m>n}

we have to cut off the

m

− n

{\displaystyle {\mathit {m}}-n}

zero padded bottom rows. and the same for the RHS: Finally, solve the square upper triangular system by back substitution:

Let's start with the solution:

import math, strformat, strutils
import arraymancer

####################################################################################################
# First part: QR decomposition.

proc eye(n: Positive): Tensor[float] =
  ## Return the (n, n) identity matrix.
  result = newTensor[float](n.int, n.int)
  for i in 0..<n: result[i, i] = 1

proc norm(v: Tensor[float]): float =
  ## return the norm of a vector.
  assert v.shape.len == 1
  result = sqrt(dot(v, v)) * sgn(v[0]).toFloat

proc houseHolder(a: Tensor[float]): Tensor[float] =
  ## return the house holder of vector "a".
  var v = a / (a[0] + norm(a))
  v[0] = 1
  result = eye(a.shape[0]) - (2 / dot(v, v)) * (v.unsqueeze(1) * v.unsqueeze(0))

proc qrDecomposition(a: Tensor): tuple[q, r: Tensor] =
  ## Return the QR decomposition of matrix "a".
  assert a.shape.len == 2
  let m = a.shape[0]
  let n = a.shape[1]
  result.q = eye(m)
  result.r = a.clone
  for i in 0..<(n - ord(m == n)):
    var h = eye(m)
    h[i..^1, i..^1] = houseHolder(result.r[i..^1, i].squeeze(1))
    result.q = result.q * h
    result.r = h * result.r

####################################################################################################
# Second part: polynomial regression example.

proc lsqr(a, b: Tensor[float]): Tensor[float] =
  let (q, r) = a.qrDecomposition()
  let n = r.shape[1]
  result = solve(r[0..<n, _], (q.transpose() * b)[0..<n])

proc polyfit(x, y: Tensor[float]; n: int): Tensor[float] =
  var z = newTensor[float](x.shape[0], n + 1)
  var t = x.reshape(x.shape[0], 1)
  for i in 0..n: z[_, i] = t^.i.toFloat
  result = lsqr(z, y.transpose())

#———————————————————————————————————————————————————————————————————————————————————————————————————

proc printMatrix(a: Tensor) =
  var str: string
  for i in 0..<a.shape[0]:
    let start = str.len
    for j in 0..<a.shape[1]:
      str.addSep(" ", start)
      str.add &"{a[i, j]:8.3f}"
    str.add '\n'
  stdout.write str

proc printVector(a: Tensor) =
  var str: string
  for i in 0..<a.shape[0]:
    str.addSep(" ")
    str.add &"{a[i]:4.1f}"
  echo str


let mat = [[12, -51,   4],
           [ 6, 167, -68],
           [-4,  24, -41]].toTensor.astype(float)

let (q, r) = mat.qrDecomposition()
echo "Q:"
printMatrix q
echo "R:"
printMatrix r
echo()

let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10].toTensor.astype(float)
let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321].toTensor.astype(float)
echo "polyfit:"
printVector polyfit(x, y, 2)