Multiple linear regression in Excel Power Query

I used BING to give me a M function. Didn't test it, but it will give you an impression on the beast you may need.
Code looks relatively clean. Actually debugging this will probably take more effort than I am able to spend ...

let
    // ------------------------------------------------------------------
    // Helper: Dot product of two lists using List.Zip
    Dot = (a as list, b as list) as number =>
        List.Sum(List.Transform(List.Zip({a, b}), each _{0} * _{1})),

    // ------------------------------------------------------------------
    // Optimized Matrix Transpose: 
    // Given a matrix (list of lists), returns its transpose.
    MatrixTranspose = (matrix as list) as list =>
        let
            numColumns = if List.IsEmpty(matrix) then 0 else List.Count(matrix{0})
        in
            List.Transform({0..numColumns - 1}, (colIndex) => List.Transform(matrix, each _{colIndex})),

    // ------------------------------------------------------------------
    // Optimized Matrix Multiply using the Dot helper and transposing B.
    MatrixMultiply = (A as list, B as list) as list =>
        let
            BT = MatrixTranspose(B)
        in
            List.Transform(A, (row) =>
                List.Transform(BT, (col) => Dot(row, col))
            ),

    // ------------------------------------------------------------------
    // Optimized Matrix Inverse using Gaussian Elimination.
    // This implementation builds an augmented matrix with the identity, 
    // then eliminates the nonpivot entries.
    MatrixInverse = (matrix as list) as list =>
        let
            n = List.Count(matrix),
            // Augment matrix with the identity matrix.
            augmented = List.Transform({0..n - 1}, (i) =>
                List.Combine({
                    matrix{i},
                    List.Transform({0..n - 1}, (j) => if i = j then 1.0 else 0.0)
                })
            ),
            // For a given pivot, normalize the pivot row and update all other rows.
            EliminatePivot = (mat as list, pivot as number) as list =>
                let
                    pivotRow = mat{pivot},
                    pivotElement = pivotRow{pivot},
                    normalizedPivotRow = List.Transform(pivotRow, each _ / pivotElement),
                    updateRow = (row as list, rowIndex as number) =>
                        if rowIndex = pivot then normalizedPivotRow
                        else
                            let
                                factor = row{pivot}
                            in
                                List.Transform({0..List.Count(row) - 1}, (j) =>
                                    row{j} - factor * normalizedPivotRow{j}
                                )
                in
                    List.Transform(
                        List.Zip({mat, List.Numbers(0, n)}),
                        each updateRow(_{0}, _{1})
                    ),
            // Eliminate each pivot sequentially.
            eliminated = List.Accumulate({0..n - 1}, augmented, (state, pivot) => EliminatePivot(state, pivot)),
            // The inverted matrix is stored in the right half of the augmented matrix.
            inverse = List.Transform(eliminated, each List.Skip(_, n))
        in
            inverse,

    // ------------------------------------------------------------------
    // MultipleRegression function computes the coefficient vector for:
    //   Y = β0 + β1 * X1 + β2 * X2 + … + ε.
    // Parameters:
    //   X – a list of lists where each inner list represents one row of predictors.
    //   Y – a list of dependent variable values.
    //   IncludeIntercept (optional) – if true (default), a column of 1's is added.
    MultipleRegression = (X as list, Y as list, Optional IncludeIntercept as nullable logical) as record =>
        let
            includeIntercept = if IncludeIntercept = null then true else IncludeIntercept,
            // If needed, prepend a column of 1's.
            Xaug = if includeIntercept then List.Transform(X, each List.InsertRange(_, 0, {1.0})) else X,
            XT = MatrixTranspose(Xaug),
            XTX = MatrixMultiply(XT, Xaug),
            XTXinv = MatrixInverse(XTX),
            // Convert Y into a column matrix.
            Ymatrix = List.Transform(Y, each { _ }),
            XTY = MatrixMultiply(XT, Ymatrix),
            Bmatrix = MatrixMultiply(XTXinv, XTY),
            // Flatten the result (each row is a one-element list).
            coefficients = List.Transform(Bmatrix, each _{0})
        in
            [Coefficients = coefficients],

    // ------------------------------------------------------------------
    // Example Usage:
    sampleX = {
        {2, 3},
        {5, 1},
        {8, 4},
        {4, 5}
    },
    sampleY = {10, 20, 30, 40},
    regressionResult = MultipleRegression(sampleX, sampleY, true)
in
    regressionResult