QR decomposition: Difference between revisions

=={{header|Ada}}==

Output matches that of Matlab solution, not tested with other matrices.

with Ada.Text_IO; use Ada.Text_IO;

with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;

Put_Line ("Q:"); Show (Q);

Put_Line ("R:"); Show (R);

{{out}}

<pre>Q:

=={{header|Axiom}}==

The following provides a generic QR decomposition for arbitrary precision floats, double floats and exact calculations:

TestPackage(R:Join(Field,RadicalCategory)): with

unitVector: NonNegativeInteger -> Vector(R)

for j in 0..n repeat

setColumn!(a,j+1,x^j)

This can be called using:

qr m

x := vector [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];

y := vector [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];

With output in exact form:

+ 6 69 58 +

[1,2,3]

The calculations are comparable to those from the default QR decomposition in R.

{{works with|BBC BASIC for Windows}}

Makes heavy use of BBC BASIC's matrix arithmetic.

@% = &2040A

INSTALL @lib$+"ARRAYLIB"

REM Create the Householder matrix H() = I - 2/vt()v() v()vt():

vvt() *= 2 / d(0) : H() = I() - vvt()

'''Output:'''

<pre>

=={{header|C}}==

#include <stdlib.h>

#include <math.h>

matrix_delete(m);

return 0;

{{out}}

<pre>

{{libheader|Math.Net}}

using MathNet.Numerics.LinearAlgebra;

using MathNet.Numerics.LinearAlgebra.Double;

Console.WriteLine(qr.R);

}

{{out}}

=={{header|C++}}==

* g++ -O3 -Wall --std=c++11 qr_standalone.cpp -o qr_standalone

*/

return EXIT_SUCCESS;

}

{{out}}

<pre>

Helper functions:

(if (zerop x)

x

C))

Main routines:

(defun make-householder (a)

(let* ((m (car (array-dimensions a)))

;; Return Q and R.

(values Q A)))

Example 1:

#2A((-0.85 0.39 0.33)

#2A((-14.0 -21.0 14.0)

( 0.0 -175.0 70.0)

Example 2, [[Polynomial regression]]:

(let* ((m (cadr (array-dimensions x)))

(A (make-array `(,m ,(+ n 1)) :initial-element 0)))

(aref x j 0))))

(aref R k k))))

(let ((x #2A((0 1 2 3 4 5 6 7 8 9 10)))

(y #2A((1 6 17 34 57 86 121 162 209 262 321))))

(polyfit x y 2))

=={{header|D}}==

{{trans|Common Lisp}}

Uses the functions copied from [[Element-wise_operations]], [[Matrix multiplication]], and [[Matrix transposition]].

std.typecons, std.numeric, std.range, std.conv;

immutable y = [[1.0, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]];

polyFit(x, y, 2).writeln;

{{out}}

<pre>[[-0.857, 0.394, 0.331],

=={{header|F_Sharp|F#}}==

// QR decomposition. Nigel Galloway: January 11th., 2022

let n=[[12.0;-51.0;4.0];[6.0;167.0;-68.0];[-4.0;24.0;-41.0]]|>MathNet.Numerics.LinearAlgebra.MatrixExtensions.matrix

let g=n|>MathNet.Numerics.LinearAlgebra.Matrix.qr

printfn $"Matrix\n------\n%A{n}\nQ\n-\n%A{g.Q}\nR\n-\n%A{g.R}"

{{out}}

<pre>

See the documentation for the [http://www.netlib.org/lapack/lapack-3.1.1/html/dgeqrf.f.html DGEQRF] and [http://www.netlib.org/lapack/lapack-3.1.1/html/dorgqr.f.html DORGQR] routines. Here the example matrix is the magic square from Albrecht Dürer's ''[https://en.wikipedia.org/wiki/Melencolia_I Melencolia I]''.

implicit none

integer, parameter :: n = 4

end do

end subroutine

{{out}}

=={{header|Futhark}}==

import "lib/github.com/diku-dk/linalg/linalg"

entry main = qr [[12.0, -51.0, 4.0],[6.0, 167.0, -68.0],[-4.0, 24.0, -41.0]]

{{out}}

<pre>

{{trans|Common Lisp}}

A fairly close port of the Common Lisp solution, this solution uses the [http://github.com/skelterjohn/go.matrix go.matrix library] for supporting functions. Note though, that go.matrix has QR decomposition, as shown in the [[Polynomial_regression#Go|Go solution]] to Polynomial regression. The solution there is coded more directly than by following the CL example here. Similarly, examination of the go.matrix QR source shows that it computes the decomposition more directly.

import (

}

return x

Output:

<pre>

===Library QR, gonum/matrix===

import (

}

return x

{{out}}

<pre>

=={{header|Haskell}}==

Square matrices only; decompose A and solve Rx = q by back substitution

import Data.List

import Text.Printf (printf)

putStrLn ("q: \n" ++ show q) >>

putStrLn ("x: \n" ++ show (backSubstitution (reverse (map reverse mR)) (reverse q) []))

{{out}}

<pre>

===QR decomposition with Numeric.LinearAlgebra===

a :: Matrix R

main = do

print $ qr a

{{out}}

<pre>((3><3)

=={{header|J}}==

h =: +@|: NB. conjugate transpose

(Q0,.Q1);(R0,.T),(-n){."1 R1

end.

+-----------------------------+----------+

| 0.857143 _0.394286 _0.331429|14 21 _14|

| 0.428571 0.902857 0.0342857| 0 175 _70|

|_0.285714 0.171429 _0.942857| 0 0 35|

'Q R'=: QR X ^/ i.3

R %.~(|:Q)+/ .* Y

'''Notes''':J offers a built-in QR decomposition function, <tt>128!:0</tt>. If J did not offer this function as a built-in, it could be written in J along the lines of the second version, which is covered in [[j:Essays/QR Decomposition|an essay on the J wiki]].

Using the [https://math.nist.gov/javanumerics/jama/ JAMA] library. Compile with: '''javac -cp Jama-1.0.3.jar Decompose.java'''.

import Jama.QRDecomposition;

qr.getR().print(10, 4);

}

{{out}}

Using the [https://dst.lbl.gov/ACSSoftware/colt/ Colt] library. Compile with: '''javac -cp colt.jar Decompose.java'''.

import cern.colt.matrix.linalg.QRDecomposition;

System.out.println(qr.getR());

}

{{out}}

Compile with: '''javac -cp commons-math3-3.6.1.jar Decompose.java'''.

import org.apache.commons.math3.linear.Array2DRowRealMatrix;

}

}

{{out}}

Using the [http://la4j.org/ la4j] library. Compile with: '''javac -cp la4j-0.6.0.jar Decompose.java'''.

import org.la4j.decomposition.QRDecompositor;

System.out.println(qr[1]);

}

{{out}}

=={{header|Julia}}==

Built-in function

{{out}}

<pre>

=={{header|Maple}}==

A:=<12,-51,4;6,167,-68;-4,24,-41>:

Q,R:=QRDecomposition(A):

Q;

{{out}}

=={{header|Mathematica}}/{{header|Wolfram Language}}==

q//MatrixForm

-> 14 21 -14

0 175 -70

=={{header|MATLAB}} / {{header|Octave}}==

6 167 -68

-4 24 -41];

Output:

<pre>Q =

=={{header|Maxima}}==

a: matrix([12, -51, 4],

4.2632564145606011E-14

load("eigen")$

unitVector(n) := ematrix(n,1,1,1,1);

a : genmatrix(lambda([i,j], if j=1 then 1.0b0 else bfloat(x[i,1]^(j-1))),

length(x),n+1),

[ 6 69 58 ]

(%o) [ 2.00000000000000000000000000002b0 ]

[ ]

=={{header|Nim}}==

{{libheader|arraymancer}}

The library “arraymancer” provides a function “qr” to get the QR decomposition. Using the Tensor type of “arraymancer” we propose here an implementation of this decomposition adapted from the Python version.

import arraymancer

let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321].toTensor.astype(float)

echo "polyfit:"

{{out}}

=={{header|PARI/GP}}==

{{works with|PARI/GP|2.6.0 and above}}

=={{header|Perl}}==

Letting the <code>PDL</code> module do all the work.

use warnings;

my ($q, $r) = mqr($a);

{{out}}

<pre>[

using matrix_mul() from [[Matrix_multiplication#Phix]]

and matrix_transpose() from [[Matrix_transposition#Phix]]

<span style="color: #000080;font-style:italic;">-- demo/rosettacode/QRdecomposition.exw</span>

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>

<span style="color: #000000;">least_squares</span><span style="color: #0000FF;">()</span>

{{out}}

<pre>

=={{header|PowerShell}}==

function qr([double[][]]$A) {

$m,$n = $A.count, $A[0].count

"X^2 X constant"

"$(polyfit $x $y 2)"

{{out}}

<pre>

{{libheader|NumPy}}

Numpy has a qr function but here is a reimplementation to show construction and use of the Householder reflections.

import numpy as np

y = np.array((1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321))

{{out}}

<pre>

=={{header|R}}==

a <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41), nrow=3, ncol=3, byrow=T)

xx <- x*x

m <- lm(y ~ x + xx)

=={{header|Racket}}==

Racket has QR-decomposition builtin:

> (require math)

> (matrix-qr (matrix [[12 -51 4]

(array #[#[6/7 -69/175 -58/175] #[3/7 158/175 6/175] #[-2/7 6/35 -33/35]])

(array #[#[14 21 -14] #[0 175 -70] #[0 0 35]])

The builtin QR-decomposition uses the Gram-Schmidt algorithm.

Here is an implementation of the Householder method:

#lang racket

(require math/matrix math/array)

[ 6 167 -68]

[-4 24 -41]]))

Output:

(array #[#[6/7 69/175 -58/175]

#[3/7 -158/175 6/175]

#[0 -175 70]

#[0 0 35]])

=={{header|Raku}}==

(formerly Perl 6)

{{Works with|rakudo|2018.06}}

use v6;

"\n",

@coef.fmt("%12.6f");

output:

This function applies the Gram Schmidt algorithm. Q is printed in the console, R can be printed or visualized.

import Prelude;

import vis::Figure;

<1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>,

<2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0>

Example using visualization

=={{header|SAS}}==

proc iml;

0 -175 70

0 0 35

=={{header|Scala}}==

{{Out}}Best seen running in your browser [https://scastie.scala-lang.org/NMueO16uQl6oivliBKZHew Scastie (remote JVM)].

import Jama.{Matrix, QRDecomposition}

print(toString(d.getR))

=={{header|SequenceL}}==

{{trans|Go}}

import <Utilities/Sequence.sl>;

import <Utilities/Conversion.sl>;

j within 1 ... n;

{{out}}

{{trans|Axiom}}

We first define a signature for a radical category joined with a field. We then define a functor with (a) structures to define operators and functions for Array and Array2, and (b) functions for the QR decomposition:

type real

val zero : real

lsqr(a,y)

end

open Real

val one = 1.0

(* output *)

=={{header|Stata}}==

See [http://www.stata.com/help.cgi?mf_qrd QR decomposition] in Stata help.

: qrd(a=(12,-51,4\6,167,-68\-4,24,-41),q=.,r=.)

2 | 0 -175 70 |

3 | 0 0 -35 |

=={{header|Tcl}}==

Assuming the presence of the Tcl solutions to these tasks: [[Element-wise operations]], [[Matrix multiplication]], [[Matrix transposition]]

{{trans|Common Lisp}}

namespace path {::tcl::mathfunc ::tcl::mathop}

proc sign x {expr {$x == 0 ? 0 : $x < 0 ? -1 : 1}}

}

return [list $Q $A]

Demonstrating:

puts "==Q=="

print_matrix [lindex $demo 0] "%f"

puts "==R=="

Output:

<pre>

=={{header|VBA}}==

Private Function vtranspose(v As Variant) As Variant

'-- transpose a vector of length m into an mx1 matrix,

Debug.Print "Q * R"

pp matrix_mul(q, r_)

<pre>A

12, -51, 4,

2, 0, 3, </pre>

Least squares

x = [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}]

y = [{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}]

Debug.Print Format(z(i), "0.#####"),

Next i

<pre>Least-squares solution: 1, 2, 3, </pre>

{{libheader|Wren-matrix}}

{{libheader|Wren-fmt}}

import "/fmt" for Fmt

Fmt.mprint(R, 8, 3)

System.print("\nQ * R:")

{{out}}

=={{header|zkl}}==

A:=GSL.Matrix(3,3).set(12.0, -51.0, 4.0,

6.0, 167.0, -68.0,

println("Q:\n",Q.format());

println("R:\n",R.format());

{{out}}

<pre>