Solve Linear Equations in Excel VBA using LU Decomposition (Faster than Matrix Inverse)
While preparing my latest post, I discovered a robust C# algorithm that solves systems of linear equations using LU Decomposition. I decided to port it to VBA and share it with you.
Why use this over MINVERSE?
Directly inverting a matrix is computationally expensive and can be prone to rounding errors. LU Decomposition factors the matrix into a Lower (L) and Upper (U) triangular matrix, which makes solving the system $Ax = B$ much faster and numerically stable.
Here is the clean, ready-to-use VBA function.
The VBA Code
Copy the code below into a standard module in your VBA Editor.
Improvements made:
-
Added Indentation: Makes the code readable and easier to debug.
-
Added Comments: Explains the three steps (Decomposition, Forward Substitution, Backward Substitution).
-
Performance: Added
Value2to read ranges into memory arrays immediately. Accessing cells inside loops is very slow; using memory arrays makes this function instant.
Option Explicit
Option Base 1
Public Function SolveLinearEquation(LeftPart As Range, RightPart As Range) As Variant
' Optimized version: Reads ranges into memory arrays for speed
Dim A As Variant, B As Variant
Dim L() As Double, U() As Double
Dim Sum As Double
Dim n As Integer, i As Integer, j As Integer, k As Integer
Dim Result() As Double, TempResult() As Double
' 1. Load Range Data into Memory Arrays
A = LeftPart.Value2
B = RightPart.Value2
' 2. Check Dimensions
n = UBound(A, 1) ' Number of rows
If UBound(A, 2) <> n Then
SolveLinearEquation = "Error: Matrix must be square"
Exit Function
End If
' 3. Initialize Arrays
ReDim L(n, n)
ReDim U(n, n)
ReDim TempResult(n)
ReDim Result(n)
' --- STEP 1: LU Decomposition ---
' Decompose Matrix A into Lower (L) and Upper (U) matrices
For i = 1 To n
' Calculate U (Upper Triangular)
For j = i To n
Sum = 0
For k = 1 To i - 1
Sum = Sum + L(i, k) * U(k, j)
Next k
U(i, j) = A(i, j) - Sum
Next j
' Calculate L (Lower Triangular)
For j = i To n
If i = j Then
L(i, i) = 1 ' Diagonal of L is 1
Else
Sum = 0
For k = 1 To i - 1
Sum = Sum + L(j, k) * U(k, i)
Next k
' Avoid division by zero
If U(i, i) = 0 Then
SolveLinearEquation = "Error: Singular Matrix"
Exit Function
End If
L(j, i) = (A(j, i) - Sum) / U(i, i)
End If
Next j
Next i
' --- STEP 2: Forward Substitution (Solve Ly = B) ---
For i = 1 To n
Sum = 0
For k = 1 To i - 1
Sum = Sum + L(i, k) * TempResult(k)
Next k
TempResult(i) = B(i, 1) - Sum
Next i
' --- STEP 3: Backward Substitution (Solve Ux = y) ---
For i = n To 1 Step -1
Sum = 0
For k = i + 1 To n
Sum = Sum + U(i, k) * Result(k)
Next k
Result(i) = (TempResult(i) - Sum) / U(i, i)
Next i
' Return the solution as a vertical array
SolveLinearEquation = Application.WorksheetFunction.Transpose(Result)
End Function
How to Use This Function
-
Open your Excel workbook and press
Alt + F11to open the VBA Editor. -
Insert a Module and paste the code above.
-
Go back to your worksheet.
Syntax: =SolveLinearEquation(Matrix_Range, Constants_Range)
Example: If your coefficient matrix is in A1:C3 (a 3×3 grid) and your solution constants are in E1:E3 (a 3×1 column), you would enter:
=SolveLinearEquation(A1:C3, E1:E3)
Important Note on Array Formulas:
-
Office 365 / Excel 2021+: Simply press Enter. The results will automatically “spill” into the cells below.
-
Older Excel Versions: Select the range of cells where you want the result (e.g., F1:F3), type the formula, and press Ctrl + Shift + Enter.
Look at this video.
Hi Rodrigue,
This matrix decomposition by triangulation comes handy for ranks over 3, although it degrades computation performance as the rank increases (it turns to be O(n³)), and some other solvers as Conjugate gradient are better suited when the rank reaches some hundreds. For information purposes only, for rank 3 (and for rank 2 also), there is a direct function to solve the inverse, as explained here https://ardoris.wordpress.com/2008/07/18/general-formula-for-the-inverse-of-a-3×3-matrix/
Kind regards mate
Here is some code Conjugate gradient from my VBA arsenal:
'CONJUGATE GRADIENT METHOD FOR SOLVING EQUATIONS
Public Sub sConjugateGradient()
Dim a() As Double
Dim B() As Double
Dim X() As Double
'-------------------------------
Dim arrA As Variant
Dim arrB As Variant
Dim row As Long
Dim col As Long
With Application
.Calculation = xlCalculationManual
.ScreenUpdating = False
End With
'arrA = Range(Cells(1, 1), Cells(xxx, xxx)).Value2
'arrB = Range(Cells(1, 1), Cells(xxx, xxx)).Value2
ReDim a(LBound(arrA, 1) To UBound(arrA, 1), LBound(arrA, 2) To UBound(arrA, 2))
ReDim B(LBound(arrB, 1) To UBound(arrB, 1))
For row = LBound(arrA, 1) To UBound(arrA, 1)
For col = LBound(arrA, 2) To UBound(arrA, 2)
a(row, col) = arrA(row, col)
Next col
B(row) = arrB(row, 1)
Next row
'-------------------------------
X() = fConjugateGradient_Solver(a(), B())
With Application
.Calculation = xlCalculationAutomatic
.ScreenUpdating = True
End With
Stop
End Sub
Public Function fConjugateGradient_Solver(ByRef a() As Double, _
ByRef B() As Double, _
Optional ByVal dbTolerance As Double = 0.000001, _
Optional ByVal iter_Max As Long = 10) As Double()
' Fletcher-Reeves algorithm (array should be symmetric)
Dim retValue As VbMsgBoxResult
Dim i As Long
Dim j As Long
Dim n As Long
Dim X() As Double
Dim g() As Double
Dim D() As Double
Dim ad() As Double
Dim iter As Long
Dim dad As Double
Dim c As Double
Dim aL As Double
Dim gg1 As Double
Dim gg2 As Double
Dim GG_Previous As Double
Dim bt As Double
'Dim dbTolerance As Double
dbTolerance = 0.000001
'Ensure matrix is square
If (UBound(a, 1) - LBound(a, 1)) _
(UBound(a, 2) - LBound(a, 2)) Then
retValue = fMsgBox("Array [A] is not square", vbCritical)
Exit Function
End If
' Solutions array
n = UBound(a, 1) - LBound(a, 1) + 1
ReDim X(LBound(a, 1) To UBound(a, 1))
' Create auxiliary arrays
ReDim g(LBound(a, 1) To UBound(a, 1))
ReDim D(LBound(a, 1) To UBound(a, 1))
ReDim ad(LBound(a, 1) To UBound(a, 1))
For i = LBound(a, 1) To UBound(a, 1)
'X(I) = 0 'Initialize array
g(i) = -B(i)
D(i) = B(i)
Next i
gg1 = 0
For i = LBound(a, 1) To UBound(a, 1)
gg1 = gg1 + (g(i) * g(i))
Next i
GG_Previous = gg1
Do While gg1 > dbTolerance
' Convergence analysis
If gg1 <= GG_Previous Then
GG_Previous = gg1
Else
' it's not converging
'GoTo ErrConvergence
End If
iter = iter + 1
dad = 0
For i = LBound(a, 1) To UBound(a, 1)
c = 0
For j = LBound(a, 1) To UBound(a, 1)
c = c + (a(i, j) * D(j))
Next j
ad(i) = c
dad = dad + (c * D(i))
Next i
'If VBA.Abs(DAD) iter_Max Then GoTo ErrSlowConvergence ' it's not converging
Loop
' Free memory
Erase g(), D(), ad() ', A(), B()
fConjugateGradient_Solver = X()
ExitProc:
Exit Function
ErrConvergence:
retValue = MsgBox("Conjugate gradient solver it's not converging, array seems to be bad conditioned" & vbNewLine & _
"[" & gg1 & ">" & GG_Previous & "] on iteration [" & iter & "]", _
vbCritical)
GoTo ExitProc
ErrSlowConvergence:
retValue = MsgBox("Conjugate gradient solver it's converging slowly, array seems to be bad conditioned" & vbNewLine & _
"Convergence = [" & gg1 & "] on iteration [" & iter & "]", _
vbCritical)
GoTo ExitProc
End Function