diff --git a/src/util.js b/src/util.js
index ecf13e87c..7b4112127 100644
--- a/src/util.js
+++ b/src/util.js
@@ -210,3 +210,258 @@ export function isInstance (arg, constructor) {
 
 	return false;
 }
+
+/**
+ * Generate a matrix of size NxN with the given diagonal values of length N.
+ *
+ * @param {number[]} values
+ * @returns {number[][]}
+ */
+export function diag (values) {
+	const n = values.length;
+	const matrix = [];
+	for (let i = 0; i < n; i++) {
+		matrix[i] = [];
+		for (let j = 0; j < n; j++) {
+			matrix[i][j] = i === j ? values[i] : 0;
+		}
+	}
+	return matrix;
+}
+
+/**
+ * Calculate the LU decomposition of an NxN matrix.
+ *
+ * P is returned as PA = UL or A = P'UL which follows Matlab and Octave opposed to Scipy which returns P as
+ * A = PUL or P'A = UL. For matrix inverse, we need P such that PA = UL and it is faster not having to invert
+ * P, even if we can invert it fairly fast as it is just a shuffled identity matrix.
+ *
+ * P is returned as a permutation matrix unless pIndices is true, in which case P would be returned as
+ * a vector containing the indexes such that A[P,:] = L*U.
+ *
+ * Reference:
+ * - https://www.statlect.com/matrix-algebra/Gaussian-elimination
+ * - https://www.sciencedirect.com/topics/mathematics/partial-pivoting
+ *
+ * @overload
+ * @param {number[][]} matrix
+ * @param {{ pIndices?: false | undefined }} [options]
+ * @returns {[number[][], number[][], number[][]]}
+ */
+/**
+ * @overload
+ * @param {number[][]} matrix
+ * @param {{ pIndices?: true }} [options]
+ * @returns {[number[], number[][], number[][]]}
+ */
+/**
+ * @param {number[][]} matrix
+ * @param {{ pIndices?: boolean | undefined }} [options]
+ * @returns {[number[] | number[][], number[][], number[][]]}
+ */
+export function lu (matrix, { pIndices = false } = {}) {
+	let p1, p2, l, u;
+
+	const n = matrix.length;
+
+	// Initialize the triangle matrices along with the permutation matrix.
+	if (pIndices) {
+		p1 = Array.from({ length: n }, (_, index) => index);
+		l = diag(new Array(n).fill(1));
+	}
+	else {
+		p2 = diag(new Array(n).fill(1));
+		l = structuredClone(p2);
+	}
+	u = structuredClone(matrix);
+
+	// Create upper and lower triangle in 'u' and 'l'. 'p' tracks the permutation (relative position of rows)
+	for (let i = 0; i < n - 1; i++) {
+		// Partial pivoting: identify the row with the maximal value in the column
+		let j = i;
+		let maximum = Math.abs(u[i][i]);
+		for (let k = i + 1; k < n; k++) {
+			const a = Math.abs(u[k][i]);
+			if (a > maximum) {
+				j = k;
+				maximum = a;
+			}
+		}
+
+		// Partial pivoting: Swap rows
+		if (j != i) {
+			// Exchange current upper triangle row with row with maximal value at pivot
+			// Update permutation matrix as well
+			[u[i], u[j]] = [u[j], u[i]];
+			if (pIndices) {
+				[p1[i], p1[j]] = [p1[j], p1[i]];
+			}
+			else {
+				[p2[i], p2[j]] = [p2[j], p2[i]];
+			}
+
+			// Only swap columns up to the pivot for the lower triangle,
+			// if on first row, there is nothing to swap
+			if (i) {
+				for (let k = 0; k < i; k++) {
+					[l[i][k], l[j][k]] = [l[j][k], l[i][k]];
+				}
+			}
+		}
+
+		// Zero at pivot point, nothing to do
+		else if (!maximum) {
+			continue;
+		}
+
+		// We have a pivot point, let's zero out everything above and below
+		// the 'l' and 'u' diagonal respectively
+		for (let j = i + 1; j < n; j++) {
+			const scalar = u[j][i] / u[i][i];
+			for (let k = i; k < n; k++) {
+				u[j][k] += -u[i][k] * scalar;
+				l[j][k] += l[i][k] * scalar;
+			}
+		}
+	}
+
+	if (pIndices) {
+		return [p1, l, u];
+	}
+	return [p2, l, u];
+}
+
+/**
+ * Forward substitution for solution of ax = b where a and b are matricies.
+ *
+ * @param {number[][]} a
+ * @param {number[][]} b
+ * @param {number} n
+ * @returns {number[][]}
+ */
+function forwardSubMatrix (a, b, n) {
+	for (let i = 0; i < n; i++) {
+		const v = b[i];
+		for (let j = 0; j < i; j++) {
+			for (let k = 0; k < n; k++) {
+				v[k] -= a[i][j] * b[j][k];
+			}
+		}
+		for (let j = 0; j < n; j++) {
+			v[j] /= a[i][i];
+		}
+	}
+	return b;
+}
+
+/**
+ * Back substitution for solution of ax = b where a and b are matricies.
+ *
+ * @param {number[][]} a
+ * @param {number[][]} b
+ * @param {number} n
+ * @returns {number[][]}
+ */
+function backSubMatrix (a, b, n) {
+	for (let i = n - 1; i > -1; i--) {
+		const v = b[i];
+		for (let j = i + 1; j < n; j++) {
+			for (var k = 0; k < n; k++) {
+				v[k] -= a[i][j] * b[j][k];
+			}
+		}
+		for (let j = 0; j < n; j++) {
+			b[i][j] /= a[i][i];
+		}
+	}
+	return b;
+}
+
+/**
+ * Forward substitution for solution of ax = b where a is a matrix and b is a vector.
+ *
+ * @param {number[][]} a
+ * @param {number[]} b
+ * @param {number} n
+ * @returns {number[]}
+ */
+function forwardSubVector (a, b, n) {
+	for (let i = 0; i < n; i++) {
+		let v = b[i];
+		for (let j = 0; j < i; j++) {
+			v -= a[i][j] * b[j];
+		}
+		b[i] = v / a[i][i];
+	}
+	return b;
+}
+
+/**
+ * Back substitution for solution of ax = b where a is a matrix and b is a vector.
+ *
+ * @param {number[][]} a
+ * @param {number[]} b
+ * @param {number} n
+ * @returns {number[]}
+ */
+function backSubVector (a, b, n) {
+	for (let i = n - 1; i > -1; i--) {
+		let v = b[i];
+		for (let j = i + 1; j < n; j++) {
+			v -= a[i][j] * b[j];
+		}
+		b[i] = v / a[i][i];
+	}
+	return b;
+}
+
+/**
+ * Invert a NxN matrix.
+ *
+ * @param {number[][]} matrix
+ * @returns {number[][]}
+ */
+export function inv (matrix) {
+	// Calculate the LU decomposition.
+	const [p, l, u] = lu(matrix);
+	const n = l.length;
+
+	// Floating point math can produce very small, non-zero determinants for singular matrices.
+	// This seems to happen in Numpy as well.
+	// Don't bother calculating sign as we only care about how close to zero we are.
+	if (l.map((row, i) => row[i] * u[i][i]).reduce((acc, val) => acc * val, 1) === 0.0) {
+		throw new Error("Matrix is singular");
+	}
+
+	// Solve for the identity matrix (will give us inverse)
+	// Permutation matrix is the identity matrix, even if shuffled.
+	return backSubMatrix(u, forwardSubMatrix(l, p, n), n);
+}
+
+/**
+ * Solve a NxN matrix representing a system of equations.
+ *
+ * @param {number[][]} matrix
+ * @param {number[]} vector
+ * @returns {number[]}
+ */
+export function solve (matrix, vector) {
+	// Calculate the LU decomposition.
+	const [p, l, u] = lu(matrix, { pIndices: true });
+	const n = l.length;
+
+	// If determinant is zero, we can't solve.
+	if (l.map((row, i) => row[i] * u[i][i]).reduce((acc, val) => acc * val, 1) === 0.0) {
+		throw new Error("Matrix is singular");
+	}
+
+	return backSubVector(
+		u,
+		forwardSubVector(
+			l,
+			p.map(i => vector[i]),
+			n,
+		),
+		n,
+	);
+}
diff --git a/test/algebra.js b/test/algebra.js
new file mode 100644
index 000000000..de6b794bb
--- /dev/null
+++ b/test/algebra.js
@@ -0,0 +1,85 @@
+import * as math from "mathjs"; // Used as test oracle
+import { solve, inv } from "../src/util.js";
+import * as check from "../node_modules/htest.dev/src/check.js";
+
+// Used to collect expected results from oracle
+function refInv (A) {
+	return math.inv(A).valueOf();
+}
+
+function refSolve (A, B) {
+	return math.lusolve(A, B).valueOf();
+}
+
+function testInvExpected (testObj) {
+	return { ...testObj, expect: refInv(...testObj.args) };
+}
+
+function testSolveExpected (testObj) {
+	return { ...testObj, expect: refSolve(...testObj.args).map(x => x[0]) };
+}
+
+const M_lin_sRGB_to_XYZ = [
+	[0.4123907992659593, 0.357584339383878, 0.1804807884018343],
+	[0.21263900587151024, 0.715168678767756, 0.07219231536073371],
+	[0.01933081871559182, 0.11919477979462598, 0.9505321522496607],
+];
+
+export default {
+	name: "Linear Algebra Tests",
+	check: check.deep(check.proximity({ epsilon: 1e-14 })),
+	tests: [
+		{
+			name: "Invert matrices",
+			tests: [
+				{
+					name: "Matrix inverse 3x3",
+					run: inv,
+					args: [M_lin_sRGB_to_XYZ],
+				},
+				{
+					name: "Matrix inverse 2x2",
+					run: inv,
+					args: [
+						[
+							[1, 2],
+							[0, 4],
+						],
+					],
+				},
+			].map(testInvExpected),
+		},
+		{
+			name: "Solve system of linear equations",
+			tests: [
+				{
+					name: "Matrix solve RGB value given an XYZ (white)",
+					run: solve,
+					args: [
+						M_lin_sRGB_to_XYZ,
+						[0.9504559270516715, 0.9999999999999999, 1.0890577507598784],
+					],
+				},
+				{
+					name: "Matrix solve RGB value given an XYZ (red)",
+					run: solve,
+					args: [
+						M_lin_sRGB_to_XYZ,
+						[0.4123907992659593, 0.21263900587151024, 0.01933081871559182],
+					],
+				},
+				{
+					name: "Solve 2D",
+					run: solve,
+					args: [
+						[
+							[1, 2],
+							[0, 4],
+						],
+						[-2, 3],
+					],
+				},
+			].map(testSolveExpected),
+		},
+	],
+};
diff --git a/test/index-fn.js b/test/index-fn.js
index 07fd33863..f0c08e2f2 100644
--- a/test/index-fn.js
+++ b/test/index-fn.js
@@ -10,6 +10,7 @@ let tests = await Promise.all(
 		"parse",
 		"contrast",
 		"multiply_matrices",
+		"algebra"
 	].map(name => import(`./${name}.js`).then(module => module.default)),
 );
 
diff --git a/test/index.json b/test/index.json
index d0421efa3..4f0e5e41d 100644
--- a/test/index.json
+++ b/test/index.json
@@ -10,5 +10,6 @@
 	"angles": "Angles",
 	"adapt": "Chromatic adaptation",
 	"multiply_matrices": "Matrix multiplication",
+	"algebra": "Linear algebra",
 	"hooks": "Hooks"
 }