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
