Mathematical formula for the OLS (least-squares) estimator:

\[\hat{\beta} = (X^T X)^{-1} X^T y\]

where one of the ways of computing the inverse of a (square) matrix is

\[A^{-1} = \frac{1}{|det(A)|} adj(A)\]

where the formula for adjoint of a matrix is given here.

However, using this method for computing the inverse of the matrix can be numerically unstable. [add cite]

Thus, a different method is used to compute the least-squares coefficients, namely through the QR-decomposition of the \(X\) matrix. This results in \(X = QR\) where Q is an orthogonal matrix and R is an upper triangular matrix. If the matrix that is used in the decomposition is \(n \times p\) where \(n > p\) (which is the case we’re concerned with here) then \[X = [Q_1 \ Q_2] \begin{bmatrix} R_1 \\ 0 \end{bmatrix}\]

where \(X = Q_1 R_1\) is known as the thin QR factorization of the data matrix X.

The inverse of a square matrix \(A\) can be written as \[A^{-1} = (QR)^{-1} = R^{-1} Q^{-1} = R^{-1} Q^T\] since \(Q\) is an orthogonal matrix. For orthogonal matrices, \[QQ^T = I \Rightarrow (Q^{-1}Q)Q^T = Q^{-1}I \Rightarrow Q^T = Q^{-1}\]

Instead of computing \((X^T X)^{-1}\) using QR and plugging it into the formula for \(\hat{\beta}\) above, we can directly use the QR matrices to compute the \(\beta\)’s using the following formula \[\hat{\beta} = R^{-1} Q^T y\]

where R^{-1} can be computed using back-substitution (R command: backsolve()) or using a different formula for computing the inverse of a triangular matrix.

The following function computes the inner product of two vectors (of equal length).

# function to compute the inner product of two vectors
# https://en.wikipedia.org/wiki/Kahan_summation_algorithm
# https://en.wikipedia.org/wiki/Dot_product#Computation
inner_prod = function(v1, v2) {
  
  stopifnot(length(v1) == length(v2))
  
  len = length(v1)
  res = vector("numeric", length = len)
  
  for (i in 1:len) {
    res[i] = v1[i] * v2[i]
  }
  
  return(sum(res))
}
inner = function(v1, v2) {
  
  stopifnot(length(v1) == length(v2))
  
  return(sum(v1 * v2))
}
# test the functions
all.equal(inner_prod(1:10, 1:10), inner(1:10, 1:10))
[1] TRUE

Using Gram-Schmidt Process

The formula/algorithm for the classical Gram-Schmidt process can be found here.

The function below implements the modified Gram-Schmidt algorithm which is a slight modification of the classical Gram-Schmidt algorithm resulting in greater numerical stability.

# Thin QR factorization and modified gram-schmidt
qr_gs = function(x) {
  
  n = nrow(x)
  p = ncol(x)
  
  q1 = matrix(0, nrow = n, ncol = p)
  r1 = matrix(0, nrow = p, ncol = p) 
  u = matrix(0, nrow = n, ncol = p)
  
  u[, 1] = x[, 1]
  
  for (k in 2:ncol(x)) {
    
    u[, k] = x[, k]
    
    # successive orthogonalization
    for (ctr in seq(1, (k - 1), 1)) {
      
      # classical GS
      #u[, k] = u[, k] - ((inner(u[, ctr], x[, k]) / inner(u[, ctr], u[, ctr])) * (u[, ctr]))
      
      # modified (stable) GS
      u[, k] = u[, k] - ((inner(u[, ctr], u[, k]) / inner(u[, ctr], u[, ctr])) * (u[, ctr]))
    }
    
  }
  
  # dividing each column by respective vector (l2) norm (rescaling each vector to have length 1)
  q1 = apply(u, 2, function(x) { x / sqrt(inner(x, x)) })
  r1 = crossprod(q1, x) # t(q1) %*% x
  
  return(list(q = q1, r = r1, u = u))
}

Check if it works: (example taken from the wikipedia article)

A = matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41), ncol = 3, byrow = TRUE)
A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41
res = qr_gs(A)
res2 = qr(A)
res$u
     [,1] [,2]  [,3]
[1,]   12  -69 -11.6
[2,]    6  158   1.2
[3,]   -4   30 -33.0
res$q
           [,1]       [,2]        [,3]
[1,]  0.8571429 -0.3942857 -0.33142857
[2,]  0.4285714  0.9028571  0.03428571
[3,] -0.2857143  0.1714286 -0.94285714
qr.Q(res2)
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714
all.equal(abs(qr.Q(res2)), abs(res$q))
[1] TRUE
round(res$r, 10)
     [,1] [,2] [,3]
[1,]   14   21  -14
[2,]    0  175  -70
[3,]    0    0   35
qr.R(res2)
     [,1] [,2] [,3]
[1,]  -14  -21   14
[2,]    0 -175   70
[3,]    0    0  -35
all.equal(abs(round(res$r, 10)), abs(qr.R(res2)))
[1] TRUE

Why are the signs from the qr call and our implementation different? The footnote on p. 20 of Wood (2006) mentions this:

In fact the QR decomposition is not uniquely defined, in that the sign of rows of Q, and corresponding columns of R, can be switched, without changing X — these sign changes are equivalent to reflections of vectors, and the sign leading to maximum numerical stability is usually selected in practice.

Check if it gives the same result for \(A^{-1}\) as R function solve()

all.equal(tcrossprod(backsolve(qr_gs(A)$r, diag(3)), qr_gs(A)$q), solve(A))
[1] TRUE

Using Householder Transformations

@TODO

Using Givens Rotations

@TODO

Computing OLS Coefficients

Taken from wikipedia.

\(\hat{\beta} = R^{-1} Q^T y\)

# respose vector y and predictor (n x p) matrix x
get_betas = function(x, y) {
  
  # get QR matrices
  res = qr_gs(x)
  
  # t(Q) %*% y
  qb = crossprod(res$q, y)
  
  # betas: R^-1 Q^T y
  coeff = as.vector(backsolve(res$r, qb))
  
  names(coeff) = colnames(x)
  
  # alternatively
  # return(as.vector(backsolve(res$r), crossprod(res$q, y)))
  
  return(coeff)
}

Test whether the beta’s are the same from the lm call and our QR implementation

# generate data
n = 100
x = matrix(rnorm(2 * n), ncol = 2)
y = 2 * x[, 1] + 3 * x[, 2] + rnorm(n)
unname(coef(lm(y ~ 0 + x))) # by lm
[1] 1.925652 2.867929
get_betas(x, y)
[1] 1.925652 2.867929

References

https://en.wikipedia.org/wiki/QR_decomposition

Simon Wood - Generalized Additive Models (2006)

Some other useful references:

LS0tCnRpdGxlOiAiQ29tcHV0aW5nIE9MUyBDb2VmZmljaWVudHMgVXNpbmcgUVIgRGVjb21wb3NpdGlvbiIKYXV0aG9yOiAiQWtzaGF0IER3aXZlZGkiCmRhdGU6ICJgciBTeXMuRGF0ZSgpYCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKTWF0aGVtYXRpY2FsIGZvcm11bGEgZm9yIHRoZSBPTFMgKGxlYXN0LXNxdWFyZXMpIGVzdGltYXRvcjoKCiQkXGhhdHtcYmV0YX0gPSAoWF5UIFgpXnstMX0gWF5UIHkkJAoKd2hlcmUgb25lIG9mIHRoZSB3YXlzIG9mIGNvbXB1dGluZyB0aGUgaW52ZXJzZSBvZiBhIChzcXVhcmUpIG1hdHJpeCBpcwoKJCRBXnstMX0gPSBcZnJhY3sxfXt8ZGV0KEEpfH0gYWRqKEEpJCQKCndoZXJlIHRoZSBmb3JtdWxhIGZvciBhZGpvaW50IG9mIGEgbWF0cml4IGlzIGdpdmVuIFtoZXJlXShodHRwczovL2VuLndpa2lwZWRpYS5vcmcvd2lraS9BZGp1Z2F0ZV9tYXRyaXgpLgoKSG93ZXZlciwgdXNpbmcgdGhpcyBtZXRob2QgZm9yIGNvbXB1dGluZyB0aGUgaW52ZXJzZSBvZiB0aGUgbWF0cml4IGNhbiBiZSBudW1lcmljYWxseSB1bnN0YWJsZS4gW2FkZCBjaXRlXQoKVGh1cywgYSBkaWZmZXJlbnQgbWV0aG9kIGlzIHVzZWQgdG8gY29tcHV0ZSB0aGUgbGVhc3Qtc3F1YXJlcyBjb2VmZmljaWVudHMsIG5hbWVseSB0aHJvdWdoIHRoZSBRUi1kZWNvbXBvc2l0aW9uIG9mIHRoZSAkWCQgbWF0cml4LiBUaGlzIHJlc3VsdHMgaW4gJFggPSBRUiQgd2hlcmUgUSBpcyBhbiBvcnRob2dvbmFsIG1hdHJpeCBhbmQgUiBpcyBhbiB1cHBlciB0cmlhbmd1bGFyIG1hdHJpeC4gSWYgdGhlIG1hdHJpeCB0aGF0IGlzIHVzZWQgaW4gdGhlIGRlY29tcG9zaXRpb24gaXMgJG4gXHRpbWVzIHAkIHdoZXJlICRuID4gcCQgKHdoaWNoIGlzIHRoZSBjYXNlIHdlJ3JlIGNvbmNlcm5lZCB3aXRoIGhlcmUpIHRoZW4gJCRYID0gW1FfMSBcIFFfMl0gXGJlZ2lue2JtYXRyaXh9IFJfMSBcXCAwIFxlbmR7Ym1hdHJpeH0kJAoKd2hlcmUgJFggPSBRXzEgUl8xJCBpcyBrbm93biBhcyB0aGUgdGhpbiBRUiBmYWN0b3JpemF0aW9uIG9mIHRoZSBkYXRhIG1hdHJpeCBYLgoKVGhlIGludmVyc2Ugb2YgYSBzcXVhcmUgbWF0cml4ICRBJCBjYW4gYmUgd3JpdHRlbiBhcyAkJEFeey0xfSA9IChRUileey0xfSA9IFJeey0xfSBRXnstMX0gPSBSXnstMX0gUV5UJCQgc2luY2UgJFEkIGlzIGFuIG9ydGhvZ29uYWwgbWF0cml4LiBGb3Igb3J0aG9nb25hbCBtYXRyaWNlcywgJCRRUV5UID0gSSBcUmlnaHRhcnJvdyAoUV57LTF9USlRXlQgPSBRXnstMX1JIFxSaWdodGFycm93IFFeVCA9IFFeey0xfSQkCgpJbnN0ZWFkIG9mIGNvbXB1dGluZyAkKFheVCBYKV57LTF9JCB1c2luZyBRUiBhbmQgcGx1Z2dpbmcgaXQgaW50byB0aGUgZm9ybXVsYSBmb3IgJFxoYXR7XGJldGF9JCBhYm92ZSwgd2UgY2FuIGRpcmVjdGx5IHVzZSB0aGUgUVIgbWF0cmljZXMgdG8gY29tcHV0ZSB0aGUgJFxiZXRhJCdzIHVzaW5nIHRoZSBmb2xsb3dpbmcgZm9ybXVsYSAkJFxoYXR7XGJldGF9ID0gUl57LTF9IFFeVCB5JCQKCndoZXJlIFJeey0xfSBjYW4gYmUgY29tcHV0ZWQgdXNpbmcgW2JhY2stc3Vic3RpdHV0aW9uXShodHRwczovL2VuLndpa2lwZWRpYS5vcmcvd2lraS9Ucmlhbmd1bGFyX21hdHJpeCNGb3J3YXJkX2FuZF9iYWNrX3N1YnN0aXR1dGlvbikgKFIgY29tbWFuZDogYGJhY2tzb2x2ZSgpYCkgb3IgdXNpbmcgW2EgZGlmZmVyZW50IGZvcm11bGEgZm9yIGNvbXB1dGluZyB0aGUgaW52ZXJzZSBvZiBhIHRyaWFuZ3VsYXIgbWF0cml4XShodHRwczovL21hdGguc3RhY2tleGNoYW5nZS5jb20vcXVlc3Rpb25zLzQ4NDEvaW52ZXJzZS1vZi1hbi1pbnZlcnRpYmxlLXRyaWFuZ3VsYXItbWF0cml4LWVpdGhlci11cHBlci1vci1sb3dlci1pcy10cmlhbmd1bGFyKS4KClRoZSBmb2xsb3dpbmcgZnVuY3Rpb24gY29tcHV0ZXMgdGhlIGlubmVyIHByb2R1Y3Qgb2YgdHdvIHZlY3RvcnMgKG9mIGVxdWFsIGxlbmd0aCkuCgpgYGB7cn0KIyBmdW5jdGlvbiB0byBjb21wdXRlIHRoZSBpbm5lciBwcm9kdWN0IG9mIHR3byB2ZWN0b3JzCiMgaHR0cHM6Ly9lbi53aWtpcGVkaWEub3JnL3dpa2kvS2FoYW5fc3VtbWF0aW9uX2FsZ29yaXRobQojIGh0dHBzOi8vZW4ud2lraXBlZGlhLm9yZy93aWtpL0RvdF9wcm9kdWN0I0NvbXB1dGF0aW9uCgppbm5lcl9wcm9kID0gZnVuY3Rpb24odjEsIHYyKSB7CiAgCiAgc3RvcGlmbm90KGxlbmd0aCh2MSkgPT0gbGVuZ3RoKHYyKSkKICAKICBsZW4gPSBsZW5ndGgodjEpCiAgcmVzID0gdmVjdG9yKCJudW1lcmljIiwgbGVuZ3RoID0gbGVuKQogIAogIGZvciAoaSBpbiAxOmxlbikgewogICAgcmVzW2ldID0gdjFbaV0gKiB2MltpXQogIH0KICAKICByZXR1cm4oc3VtKHJlcykpCn0KCmlubmVyID0gZnVuY3Rpb24odjEsIHYyKSB7CiAgCiAgc3RvcGlmbm90KGxlbmd0aCh2MSkgPT0gbGVuZ3RoKHYyKSkKICAKICByZXR1cm4oc3VtKHYxICogdjIpKQp9CgojIHRlc3QgdGhlIGZ1bmN0aW9ucwphbGwuZXF1YWwoaW5uZXJfcHJvZCgxOjEwLCAxOjEwKSwgaW5uZXIoMToxMCwgMToxMCkpCmBgYAoKIyMjIFVzaW5nIEdyYW0tU2NobWlkdCBQcm9jZXNzCgpUaGUgZm9ybXVsYS9hbGdvcml0aG0gZm9yIHRoZSBjbGFzc2ljYWwgR3JhbS1TY2htaWR0IHByb2Nlc3MgY2FuIGJlIGZvdW5kIFtoZXJlXShodHRwczovL2VuLndpa2lwZWRpYS5vcmcvd2lraS9RUl9kZWNvbXBvc2l0aW9uI1VzaW5nX3RoZV9HcmFtLkUyLjgwLjkzU2NobWlkdF9wcm9jZXNzKS4KClRoZSBmdW5jdGlvbiBiZWxvdyBpbXBsZW1lbnRzIHRoZSBbbW9kaWZpZWQgR3JhbS1TY2htaWR0IGFsZ29yaXRobV0oaHR0cHM6Ly9lbi53aWtpcGVkaWEub3JnL3dpa2kvR3JhbSVFMiU4MCU5M1NjaG1pZHRfcHJvY2VzcyNOdW1lcmljYWxfc3RhYmlsaXR5KSB3aGljaCBpcyBhIHNsaWdodCBtb2RpZmljYXRpb24gb2YgdGhlIGNsYXNzaWNhbCBHcmFtLVNjaG1pZHQgYWxnb3JpdGhtIHJlc3VsdGluZyBpbiBncmVhdGVyIG51bWVyaWNhbCBzdGFiaWxpdHkuCgpgYGB7cn0KIyBUaGluIFFSIGZhY3Rvcml6YXRpb24gYW5kIG1vZGlmaWVkIGdyYW0tc2NobWlkdApxcl9ncyA9IGZ1bmN0aW9uKHgpIHsKICAKICBuID0gbnJvdyh4KQogIHAgPSBuY29sKHgpCiAgCiAgcTEgPSBtYXRyaXgoMCwgbnJvdyA9IG4sIG5jb2wgPSBwKQogIHIxID0gbWF0cml4KDAsIG5yb3cgPSBwLCBuY29sID0gcCkgCiAgdSA9IG1hdHJpeCgwLCBucm93ID0gbiwgbmNvbCA9IHApCiAgCiAgdVssIDFdID0geFssIDFdCiAgCiAgZm9yIChrIGluIDI6bmNvbCh4KSkgewogICAgCiAgICB1Wywga10gPSB4Wywga10KICAgIAogICAgIyBzdWNjZXNzaXZlIG9ydGhvZ29uYWxpemF0aW9uCiAgICBmb3IgKGN0ciBpbiBzZXEoMSwgKGsgLSAxKSwgMSkpIHsKICAgICAgCiAgICAgICMgY2xhc3NpY2FsIEdTCiAgICAgICN1Wywga10gPSB1Wywga10gLSAoKGlubmVyKHVbLCBjdHJdLCB4Wywga10pIC8gaW5uZXIodVssIGN0cl0sIHVbLCBjdHJdKSkgKiAodVssIGN0cl0pKQogICAgICAKICAgICAgIyBtb2RpZmllZCAoc3RhYmxlKSBHUwogICAgICB1Wywga10gPSB1Wywga10gLSAoKGlubmVyKHVbLCBjdHJdLCB1Wywga10pIC8gaW5uZXIodVssIGN0cl0sIHVbLCBjdHJdKSkgKiAodVssIGN0cl0pKQogICAgfQogICAgCiAgfQogIAogICMgZGl2aWRpbmcgZWFjaCBjb2x1bW4gYnkgcmVzcGVjdGl2ZSB2ZWN0b3IgKGwyKSBub3JtIChyZXNjYWxpbmcgZWFjaCB2ZWN0b3IgdG8gaGF2ZSBsZW5ndGggMSkKICBxMSA9IGFwcGx5KHUsIDIsIGZ1bmN0aW9uKHgpIHsgeCAvIHNxcnQoaW5uZXIoeCwgeCkpIH0pCiAgcjEgPSBjcm9zc3Byb2QocTEsIHgpICMgdChxMSkgJSolIHgKICAKICByZXR1cm4obGlzdChxID0gcTEsIHIgPSByMSwgdSA9IHUpKQp9CmBgYAoKQ2hlY2sgaWYgaXQgd29ya3M6IChleGFtcGxlIHRha2VuIGZyb20gdGhlIHdpa2lwZWRpYSBhcnRpY2xlKQoKYGBge3J9CkEgPSBtYXRyaXgoYygxMiwgLTUxLCA0LCA2LCAxNjcsIC02OCwgLTQsIDI0LCAtNDEpLCBuY29sID0gMywgYnlyb3cgPSBUUlVFKQpBCgpyZXMgPSBxcl9ncyhBKQpyZXMyID0gcXIoQSkKCnJlcyR1CnJlcyRxCnFyLlEocmVzMikKYWxsLmVxdWFsKGFicyhxci5RKHJlczIpKSwgYWJzKHJlcyRxKSkKcm91bmQocmVzJHIsIDEwKQpxci5SKHJlczIpCmFsbC5lcXVhbChhYnMocm91bmQocmVzJHIsIDEwKSksIGFicyhxci5SKHJlczIpKSkKYGBgCgpXaHkgYXJlIHRoZSBzaWducyBmcm9tIHRoZSBgcXJgIGNhbGwgYW5kIG91ciBpbXBsZW1lbnRhdGlvbiBkaWZmZXJlbnQ/IFRoZSBmb290bm90ZSBvbiBwLiAyMCBvZiBXb29kICgyMDA2KSBtZW50aW9ucyB0aGlzOgoKYGBgCkluIGZhY3QgdGhlIFFSIGRlY29tcG9zaXRpb24gaXMgbm90IHVuaXF1ZWx5IGRlZmluZWQsIGluIHRoYXQgdGhlIHNpZ24gb2Ygcm93cyBvZiBRLCBhbmQgY29ycmVzcG9uZGluZyBjb2x1bW5zIG9mIFIsIGNhbiBiZSBzd2l0Y2hlZCwgd2l0aG91dCBjaGFuZ2luZyBYIOKAlCB0aGVzZSBzaWduIGNoYW5nZXMgYXJlIGVxdWl2YWxlbnQgdG8gcmVmbGVjdGlvbnMgb2YgdmVjdG9ycywgYW5kIHRoZSBzaWduIGxlYWRpbmcgdG8gbWF4aW11bSBudW1lcmljYWwgc3RhYmlsaXR5IGlzIHVzdWFsbHkgc2VsZWN0ZWQgaW4gcHJhY3RpY2UuCmBgYApDaGVjayBpZiBpdCBnaXZlcyB0aGUgc2FtZSByZXN1bHQgZm9yICRBXnstMX0kIGFzIFIgZnVuY3Rpb24gYHNvbHZlKClgCgpgYGB7cn0KYWxsLmVxdWFsKHRjcm9zc3Byb2QoYmFja3NvbHZlKHFyX2dzKEEpJHIsIGRpYWcoMykpLCBxcl9ncyhBKSRxKSwgc29sdmUoQSkpCmBgYAoKIyMjIFVzaW5nIEhvdXNlaG9sZGVyIFRyYW5zZm9ybWF0aW9ucwoKQFRPRE8KCiMjIyBVc2luZyBHaXZlbnMgUm90YXRpb25zCgpAVE9ETwoKIyMjIENvbXB1dGluZyBPTFMgQ29lZmZpY2llbnRzCgpUYWtlbiBmcm9tIFt3aWtpcGVkaWFdKGh0dHBzOi8vZW4ud2lraXBlZGlhLm9yZy93aWtpL1FSX2RlY29tcG9zaXRpb24jVXNpbmdfZm9yX3NvbHV0aW9uX3RvX2xpbmVhcl9pbnZlcnNlX3Byb2JsZW1zKS4KCiRcaGF0e1xiZXRhfSA9IFJeey0xfSBRXlQgeSQKCmBgYHtyfQojIHJlc3Bvc2UgdmVjdG9yIHkgYW5kIHByZWRpY3RvciAobiB4IHApIG1hdHJpeCB4CmdldF9iZXRhcyA9IGZ1bmN0aW9uKHgsIHkpIHsKICAKICAjIGdldCBRUiBtYXRyaWNlcwogIHJlcyA9IHFyX2dzKHgpCiAgCiAgIyB0KFEpICUqJSB5CiAgcWIgPSBjcm9zc3Byb2QocmVzJHEsIHkpCiAgCiAgIyBiZXRhczogUl4tMSBRXlQgeQogIGNvZWZmID0gYXMudmVjdG9yKGJhY2tzb2x2ZShyZXMkciwgcWIpKQogIAogIG5hbWVzKGNvZWZmKSA9IGNvbG5hbWVzKHgpCiAgCiAgIyBhbHRlcm5hdGl2ZWx5CiAgIyByZXR1cm4oYXMudmVjdG9yKGJhY2tzb2x2ZShyZXMkciksIGNyb3NzcHJvZChyZXMkcSwgeSkpKQogIAogIHJldHVybihjb2VmZikKfQpgYGAKClRlc3Qgd2hldGhlciB0aGUgYmV0YSdzIGFyZSB0aGUgc2FtZSBmcm9tIHRoZSBgbG1gIGNhbGwgYW5kIG91ciBRUiBpbXBsZW1lbnRhdGlvbgoKYGBge3J9CiMgZ2VuZXJhdGUgZGF0YQpuID0gMTAwCnggPSBtYXRyaXgocm5vcm0oMiAqIG4pLCBuY29sID0gMikKeSA9IDIgKiB4WywgMV0gKyAzICogeFssIDJdICsgcm5vcm0obikKCnVubmFtZShjb2VmKGxtKHkgfiAwICsgeCkpKSAjIGJ5IGxtCgpnZXRfYmV0YXMoeCwgeSkKYGBgCgpgYGB7ciwgaW5jbHVkZT1GQUxTRSwgZXZhbD1GQUxTRX0Kb3V0ID0gbWF0cml4KDAsIG5yb3cgPSBucm93KEEpLCBuY29sID0gbmNvbChBKSkKZm9yIChpIGluIDE6bmNvbChBKSkgewogIGZvciAoaiBpbiAxOm5jb2woQSkpIHsKICAgIG91dFtpLCBqXSA9IHN1bShBW2ksIF0gKiBBWywgal0pCiAgfQp9Cm91dApBICUqJSBBCmFsbC5lcXVhbChvdXQsIEEgJSolIEEpCmBgYAoKIyMjIFJlZmVyZW5jZXMKCmh0dHBzOi8vZW4ud2lraXBlZGlhLm9yZy93aWtpL1FSX2RlY29tcG9zaXRpb24KClNpbW9uIFdvb2QgLSBHZW5lcmFsaXplZCBBZGRpdGl2ZSBNb2RlbHMgKDIwMDYpCgpTb21lIG90aGVyIHVzZWZ1bCByZWZlcmVuY2VzOgoKKiBodHRwczovL3N0YWNrb3ZlcmZsb3cuY29tL3F1ZXN0aW9ucy8zOTg0OTk0MS93cml0aW5nLWEtaG91c2Vob2xkZXItcXItZmFjdG9yaXphdGlvbi1mdW5jdGlvbi1pbi1yLWNvZGUKCiogaHR0cHM6Ly9zdGFja292ZXJmbG93LmNvbS9xdWVzdGlvbnMvMzg3NDExMjcvbG0td2hhdC1pcy1xcmF1eC1yZXR1cm5lZC1ieS1xci1kZWNvbXBvc2l0aW9uLWluLWxpbnBhY2stbGFwYWNrCgoqIGh0dHBzOi8vc3RhY2tvdmVyZmxvdy5jb20vcXVlc3Rpb25zLzEzMjQ4OTYwL3FyLWRlY29tcG9zaXRpb24tZGlmZmVyZW50LWluLWxtLWFuZC1iaWdsbQo=