diff --git a/src/FSharp.Stats/FSharp.Stats.fsproj b/src/FSharp.Stats/FSharp.Stats.fsproj
index 633a4162..c7fc980a 100644
--- a/src/FSharp.Stats/FSharp.Stats.fsproj
+++ b/src/FSharp.Stats/FSharp.Stats.fsproj
@@ -59,6 +59,7 @@
     <Compile Include="SpecialFunctions\Logistic.fs" />
     <Compile Include="SpecialFunctions\Binomial.fs" />
     <!-- Algebra -->
+    <Compile Include="LinearAlgebra.fs" />
     <!-- RootFinding -->
     <Compile Include="RootFinding\Brent.fs" />
     <!-- Integration -->
diff --git a/src/FSharp.Stats/LinearAlgebra.fs b/src/FSharp.Stats/LinearAlgebra.fs
new file mode 100644
index 00000000..2ef0b374
--- /dev/null
+++ b/src/FSharp.Stats/LinearAlgebra.fs
@@ -0,0 +1,90 @@
+namespace FSharp.Stats
+
+open FsMath
+open FsMath.Algebra
+
+/// <summary>
+///   Provides matrix decomposition algorithms for linear algebra operations.
+/// </summary>
+module LinearAlgebra =
+
+    /// <summary>
+    ///   Specifies the algorithm used for QR decomposition.
+    /// </summary>
+    /// <remarks>
+    ///   The two algorithms differ in the shape of the output matrices:
+    ///   <list type="bullet">
+    ///     <item><description>
+    ///       <b>Householder</b>: full (uneconomised) decomposition.
+    ///       For an m×n input, produces Q (m×m) and R (m×n).
+    ///     </description></item>
+    ///     <item><description>
+    ///       <b>GramSchmidt</b>: thin (economy) decomposition.
+    ///       For an m×n input with m ≥ n, produces Q (m×n) and R (n×n).
+    ///       Useful when only the column space of A is needed and memory is
+    ///       important, or when you want Q to be exactly m×n rather than m×m.
+    ///     </description></item>
+    ///   </list>
+    ///   Both satisfy A = Q * R; both produce an orthonormal Q (Q^T Q = I).
+    /// </remarks>
+    type QRMethod =
+        /// Full decomposition via Householder reflections. Q is m×m, R is m×n.
+        | Householder
+        /// Thin (economy) decomposition via modified Gram-Schmidt orthogonalisation. Q is m×n, R is n×n.
+        | GramSchmidt
+
+    /// <summary>
+    ///   QR decomposition of a matrix A into an orthogonal matrix Q and an upper-triangular matrix R such that A = Q * R.
+    /// </summary>
+    module QR =
+
+        /// <summary>
+        ///   Decomposes matrix A into Q * R using the specified method.
+        /// </summary>
+        /// <param name="method">
+        ///   <see cref="QRMethod.Householder"/> for the full decomposition (Q is m×m);
+        ///   <see cref="QRMethod.GramSchmidt"/> for the thin/economy decomposition (Q is m×n).
+        /// </param>
+        /// <param name="A">The input matrix to decompose. Must have at least as many rows as columns for Gram-Schmidt.</param>
+        /// <returns>A tuple (Q, R) satisfying A = Q * R.</returns>
+        /// <example>
+        /// <code>
+        ///   open FsMath                       // for matrix / vector literals
+        ///   open FSharp.Stats.LinearAlgebra
+        ///
+        ///   let A = matrix [[12.;-51.;4.];[6.;167.;-68.];[-4.;24.;-41.]]
+        ///
+        ///   // Gram-Schmidt – thin Q (3×3 for a square input)
+        ///   let (qGS, rGS) = QR.decompose GramSchmidt A
+        ///
+        ///   // Householder – full Q (3×3 for a square input)
+        ///   let (qHH, rHH) = QR.decompose Householder A
+        ///
+        ///   // For a 4×3 input the difference is more visible:
+        ///   let B = matrix [[1.;2.];[3.;4.];[5.;6.];[7.;8.]]
+        ///   let (qGS4x2, rGS2x2) = QR.decompose GramSchmidt B   // Q: 4×2, R: 2×2
+        ///   let (qHH4x4, rHH4x2) = QR.decompose Householder B   // Q: 4×4, R: 4×2
+        /// </code>
+        /// </example>
+        let decompose (method: QRMethod) (A: Matrix<float>) : Matrix<float> * Matrix<float> =
+            match method with
+            | Householder -> LinearAlgebra.qrDecompose A
+            | GramSchmidt -> LinearAlgebra.qrModifiedGramSchmidt A
+
+        /// <summary>
+        ///   Decomposes matrix A using Householder reflections (full QR).
+        ///   For an m×n matrix, Q is m×m and R is m×n.
+        /// </summary>
+        /// <param name="A">The input matrix to decompose.</param>
+        /// <returns>A tuple (Q, R) satisfying A = Q * R.</returns>
+        let householder (A: Matrix<float>) : Matrix<float> * Matrix<float> =
+            LinearAlgebra.qrDecompose A
+
+        /// <summary>
+        ///   Decomposes matrix A using modified Gram-Schmidt orthogonalisation (thin/economy QR).
+        ///   For an m×n matrix with m ≥ n, Q is m×n and R is n×n.
+        /// </summary>
+        /// <param name="A">The input matrix to decompose. m must be ≥ n.</param>
+        /// <returns>A tuple (Q, R) satisfying A = Q * R.</returns>
+        let gramSchmidt (A: Matrix<float>) : Matrix<float> * Matrix<float> =
+            LinearAlgebra.qrModifiedGramSchmidt A
diff --git a/tests/FSharp.Stats.Tests/Fitting.fs b/tests/FSharp.Stats.Tests/Fitting.fs
index 91111c40..a4d08a8a 100644
--- a/tests/FSharp.Stats.Tests/Fitting.fs
+++ b/tests/FSharp.Stats.Tests/Fitting.fs
@@ -55,6 +55,95 @@ let leastSquaresCholeskyTests =
         )
     ]
 open FSharp.Stats.Fitting.Spline
+[<Tests>]
+let qrDecompositionTests =
+    // Classic Golub-Van Loan 3×3 example: known Q and R
+    let A = matrix [[12.;-51.;4.];[6.;167.;-68.];[-4.;24.;-41.]]
+    // Rectangular tall matrix for thin-QR tests
+    let B = matrix [[1.;2.];[3.;4.];[5.;6.];[7.;8.]]
+
+    let matClose (m1: Matrix<float>) (m2: Matrix<float>) label =
+        Expect.equal m1.NumRows m2.NumRows (label + " – row count")
+        Expect.equal m1.NumCols m2.NumCols (label + " – col count")
+        for i in 0 .. m1.NumRows - 1 do
+            for j in 0 .. m1.NumCols - 1 do
+                Expect.floatClose Accuracy.medium m1.[i,j] m2.[i,j] (sprintf "%s [%d,%d]" label i j)
+
+    testList "QR Decomposition" [
+
+        testCase "GramSchmidt: A = Q * R (square)" (fun () ->
+            let (q, r) = LinearAlgebra.QR.gramSchmidt A
+            let qr = q * r
+            matClose qr A "GS Q*R should equal A"
+        )
+
+        testCase "Householder: A = Q * R (square)" (fun () ->
+            let (q, r) = LinearAlgebra.QR.householder A
+            let qr = q * r
+            matClose qr A "Householder Q*R should equal A"
+        )
+
+        testCase "GramSchmidt: Q is orthonormal (Q^T Q = I)" (fun () ->
+            let (q, _) = LinearAlgebra.QR.gramSchmidt A
+            let qtq = q.Transpose() * q
+            let eye = Matrix.identity q.NumCols
+            matClose qtq eye "GS Q^T Q should be identity"
+        )
+
+        testCase "Householder: Q is orthogonal (Q^T Q = I)" (fun () ->
+            let (q, _) = LinearAlgebra.QR.householder A
+            let qtq = q.Transpose() * q
+            let eye = Matrix.identity q.NumCols
+            matClose qtq eye "Householder Q^T Q should be identity"
+        )
+
+        testCase "GramSchmidt: R is upper triangular (square)" (fun () ->
+            let (_, r) = LinearAlgebra.QR.gramSchmidt A
+            for i in 1 .. r.NumRows - 1 do
+                for j in 0 .. i - 1 do
+                    Expect.floatClose Accuracy.medium r.[i,j] 0. (sprintf "GS R[%d,%d] should be 0" i j)
+        )
+
+        testCase "Householder: R is upper triangular (square)" (fun () ->
+            let (_, r) = LinearAlgebra.QR.householder A
+            for i in 1 .. r.NumRows - 1 do
+                for j in 0 .. i - 1 do
+                    Expect.floatClose Accuracy.medium r.[i,j] 0. (sprintf "HH R[%d,%d] should be 0" i j)
+        )
+
+        testCase "GramSchmidt thin QR: dimensions for 4×2 input" (fun () ->
+            let (q, r) = LinearAlgebra.QR.gramSchmidt B
+            Expect.equal q.NumRows 4 "GS thin Q should have 4 rows"
+            Expect.equal q.NumCols 2 "GS thin Q should have 2 cols"
+            Expect.equal r.NumRows 2 "GS thin R should have 2 rows"
+            Expect.equal r.NumCols 2 "GS thin R should have 2 cols"
+        )
+
+        testCase "Householder full QR: dimensions for 4×2 input" (fun () ->
+            let (q, r) = LinearAlgebra.QR.householder B
+            Expect.equal q.NumRows 4 "HH full Q should have 4 rows"
+            Expect.equal q.NumCols 4 "HH full Q should have 4 cols"
+            Expect.equal r.NumRows 4 "HH full R should have 4 rows"
+            Expect.equal r.NumCols 2 "HH full R should have 2 cols"
+        )
+
+        testCase "GramSchmidt thin QR: A = Q * R (rectangular)" (fun () ->
+            let (q, r) = LinearAlgebra.QR.gramSchmidt B
+            let qr = q * r
+            matClose qr B "GS thin Q*R should equal B"
+        )
+
+        testCase "decompose dispatches correctly" (fun () ->
+            let (qGS, _) = LinearAlgebra.QR.decompose LinearAlgebra.GramSchmidt A
+            let (qGS2, _) = LinearAlgebra.QR.gramSchmidt A
+            matClose qGS (qGS2 |> id) "decompose GramSchmidt should match gramSchmidt"
+            let (qHH, _) = LinearAlgebra.QR.decompose LinearAlgebra.Householder A
+            let (qHH2, _) = LinearAlgebra.QR.householder A
+            matClose qHH qHH2 "decompose Householder should match householder"
+        )
+    ]
+
+
 [<Tests>]
 let splineTests = 
     testList "Fitting.Spline" [