diff --git a/README.md b/README.md
index 3426413..985123c 100644
--- a/README.md
+++ b/README.md
@@ -11,10 +11,13 @@ A linear algebra library written in TypeScript that focuses on generality, exten
 - Matrix exponentials
 - Elementary row operations
 - Gauss-Jordan elimination
+- Matrix Factorizations:
+  - Cholesky Decomposition
+  - LU Decomposition
+  - QR Decomposition
+  - Singular Value Decomposition
 - Differentiation via finite differences
 - Least-Squares Regression for arbitrary model functions
-- LU Decomposition
-- QR Decomposition
 - Eigenvalue / Eigenvector finding
 - And more to come!
 
diff --git a/benchmarks/benchmark.ts b/benchmarks/benchmark.ts
index 3fc0800..9b08f32 100644
--- a/benchmarks/benchmark.ts
+++ b/benchmarks/benchmark.ts
@@ -1,7 +1,8 @@
 /* tslint:disable:no-console */
 import * as fs from 'fs';
 import {
-  calculateLinearLeastSquaresApproximation,
+  calculateLinearLeastSquares,
+  calculateSingularValueDecomposition,
   inverse,
   linspace,
   NumberMatrix,
@@ -9,7 +10,7 @@ import {
   solveByGaussianElimination
 } from '../src';
 
-const repetitions: number = 5;
+const repetitions = 5;
 
 interface Benchmark {
   description: string;
@@ -25,7 +26,7 @@ interface BenchmarkEntry {
 
 generateBenchmarks();
 
-function generateBenchmarks() {
+function generateBenchmarks(): void {
   generateBenchmark(
     'matrix-plus-matrix',
     1,
@@ -91,7 +92,17 @@ function generateBenchmarks() {
         .fromColumnVectors([xs, ys])
         .getRowVectors();
     },
-    dataPoints => calculateLinearLeastSquaresApproximation(dataPoints)
+    dataPoints => calculateLinearLeastSquares(dataPoints)
+  );
+
+  generateBenchmark(
+    'singular-value-decomposition',
+    1,
+    32,
+    N => ({
+      A: NumberMatrix.builder().random(N)
+    }),
+    ({ A }) => calculateSingularValueDecomposition(A)
   );
 }
 
@@ -101,7 +112,7 @@ function generateBenchmark<T>(
   maxSize: number,
   prepare: (N: number) => T,
   runProcess: (data: T) => void
-) {
+): void {
   console.log(`Benchmarking ${description}...`);
 
   const benchmark: Benchmark = {
@@ -127,7 +138,7 @@ function generateBenchmark<T>(
   writeBenchmarkToFile(benchmark);
 }
 
-function generateEntry(size: number, timings: number[]) {
+function generateEntry(size: number, timings: number[]): BenchmarkEntry {
   let minTime = Number.MAX_VALUE;
   let maxTime = -Number.MAX_VALUE;
   let total = 0;
@@ -147,7 +158,7 @@ function generateEntry(size: number, timings: number[]) {
   return { size, minTime, maxTime, meanTime };
 }
 
-function writeBenchmarkToFile(benchmark: Benchmark) {
+function writeBenchmarkToFile(benchmark: Benchmark): void {
   const csvHeader = `Size,Min Time,Max Time,Mean Time\n`;
   const csvBody = benchmark.entries
     .map(entry => `${entry.size},${entry.minTime},${entry.maxTime},${entry.meanTime}`)
diff --git a/benchmarks/latest/singular-value-decomposition.csv b/benchmarks/latest/singular-value-decomposition.csv
new file mode 100644
index 0000000..a02f2b9
--- /dev/null
+++ b/benchmarks/latest/singular-value-decomposition.csv
@@ -0,0 +1,33 @@
+Size,Min Time,Max Time,Mean Time
+1,0,3,1.2
+2,2,5,3.2
+3,2,5,3.6
+4,4,8,5.4
+5,7,9,8.2
+6,11,12,11.8
+7,15,20,17.2
+8,23,26,24.2
+9,32,39,34
+10,45,46,45.6
+11,59,63,60.6
+12,81,88,85
+13,104,106,105
+14,132,137,134
+15,165,179,169.4
+16,205,219,211.2
+17,251,266,255.6
+18,311,325,316
+19,380,394,384.6
+20,460,469,464.8
+21,542,554,547.6
+22,682,696,688.8
+23,800,834,810.6
+24,935,990,951.4
+25,1092,1105,1099.2
+26,1269,1287,1277.4
+27,1464,1482,1471.2
+28,1679,1694,1685.4
+29,1962,2229,2113
+30,2233,2843,2520.2
+31,2482,2793,2586.4
+32,2807,3444,3173.4
\ No newline at end of file
diff --git a/docs/VECTOR.api.md b/docs/VECTOR.api.md
index 76b8041..a314532 100644
--- a/docs/VECTOR.api.md
+++ b/docs/VECTOR.api.md
@@ -67,6 +67,9 @@ export abstract class ArrayVector<S> implements Vector<S> {
 // @public
 export function backwardDifferenceMatrix(binCount: number): NumberMatrix;
 
+// @public
+export function calculateCholeskyDecomposition<S>(A: Matrix<S>): CholeskyDecomposition<S> | undefined;
+
 // @public
 export function calculateEigenvalues<S>(A: Matrix<S>, numIterations?: number, throwOnFailure?: boolean): Vector<S>;
 
@@ -82,9 +85,18 @@ export function calculateLUDecomposition<S>(A: Matrix<S>): LUDecomposition<S>;
 // @public
 export function calculateQRDecomposition<S>(A: Matrix<S>): QRDecomposition<S>;
 
+// @public
+export function calculateSingularValueDecomposition<S>(A: Matrix<S>): SingularValueDecomposition<S>;
+
 // @public
 export function centralDifferenceMatrix(binCount: number): NumberMatrix;
 
+// @public
+export interface CholeskyDecomposition<S> {
+    // (undocumented)
+    L: Matrix<S>;
+}
+
 // @public
 export function columnSumSupremumNorm<S>(A: Matrix<S>): number;
 
@@ -425,6 +437,9 @@ export interface QRDecomposition<S> {
     R: Matrix<S>;
 }
 
+// @public
+export function rank<S>(matrix: Matrix<S>): number;
+
 // @public
 export function reducedRowEchelonForm<S>(matrix: Matrix<S>): Matrix<S>;
 
@@ -465,6 +480,16 @@ export abstract class ScalarOperations<S> {
     zero(): S;
 }
 
+// @public
+export interface SingularValueDecomposition<S> {
+    // (undocumented)
+    Sigma: Matrix<S>;
+    // (undocumented)
+    U: Matrix<S>;
+    // (undocumented)
+    V: Matrix<S>;
+}
+
 // Warning: (ae-forgotten-export) The symbol "LinearSolution" needs to be exported by the entry point index.d.ts
 // 
 // @public
diff --git a/docs/api/vector.calculatecholeskydecomposition.md b/docs/api/vector.calculatecholeskydecomposition.md
new file mode 100644
index 0000000..1042994
--- /dev/null
+++ b/docs/api/vector.calculatecholeskydecomposition.md
@@ -0,0 +1,30 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [calculateCholeskyDecomposition](./vector.calculatecholeskydecomposition.md)
+
+## calculateCholeskyDecomposition() function
+
+Uses the serial version of the Cholesky algorith to calculate the Cholesky decomposition of a matrix `A`<!-- -->.
+
+<b>Signature:</b>
+
+```typescript
+export declare function calculateCholeskyDecomposition<S>(A: Matrix<S>): CholeskyDecomposition<S> | undefined;
+```
+
+## Parameters
+
+|  Parameter | Type | Description |
+|  --- | --- | --- |
+|  A | <code>Matrix&lt;S&gt;</code> | The matrix to decompose |
+
+<b>Returns:</b>
+
+`CholeskyDecomposition<S> | undefined`
+
+## Remarks
+
+A Cholesky decomposition of a matrix `A` consists of a lower-triangular matrix `L` such that \_LL\* = A\_.
+
+A Cholesky decomposition only exists if `A` is symmetrix and positive-definite.
+
diff --git a/docs/api/vector.calculatesingularvaluedecomposition.md b/docs/api/vector.calculatesingularvaluedecomposition.md
new file mode 100644
index 0000000..1b2e176
--- /dev/null
+++ b/docs/api/vector.calculatesingularvaluedecomposition.md
@@ -0,0 +1,28 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [calculateSingularValueDecomposition](./vector.calculatesingularvaluedecomposition.md)
+
+## calculateSingularValueDecomposition() function
+
+Uses the Power Method to calculate the Singular Value Decomposition of a matrix `A`
+
+<b>Signature:</b>
+
+```typescript
+export declare function calculateSingularValueDecomposition<S>(A: Matrix<S>): SingularValueDecomposition<S>;
+```
+
+## Parameters
+
+|  Parameter | Type | Description |
+|  --- | --- | --- |
+|  A | <code>Matrix&lt;S&gt;</code> | the matrix to decompose |
+
+<b>Returns:</b>
+
+`SingularValueDecomposition<S>`
+
+## Remarks
+
+A Singular Value Decomposition consists of orthogonal matrices `U` and `V` and a diagonal matrix `Sigma` such that \_USigmaV\* = A\_
+
diff --git a/docs/api/vector.choleskydecomposition.l.md b/docs/api/vector.choleskydecomposition.l.md
new file mode 100644
index 0000000..1a8470c
--- /dev/null
+++ b/docs/api/vector.choleskydecomposition.l.md
@@ -0,0 +1,11 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [CholeskyDecomposition](./vector.choleskydecomposition.md) &gt; [L](./vector.choleskydecomposition.l.md)
+
+## CholeskyDecomposition.L property
+
+<b>Signature:</b>
+
+```typescript
+L: Matrix<S>;
+```
diff --git a/docs/api/vector.choleskydecomposition.md b/docs/api/vector.choleskydecomposition.md
new file mode 100644
index 0000000..328621b
--- /dev/null
+++ b/docs/api/vector.choleskydecomposition.md
@@ -0,0 +1,20 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [CholeskyDecomposition](./vector.choleskydecomposition.md)
+
+## CholeskyDecomposition interface
+
+The result of a Cholesky Decomposition
+
+<b>Signature:</b>
+
+```typescript
+export interface CholeskyDecomposition<S> 
+```
+
+## Properties
+
+|  Property | Type | Description |
+|  --- | --- | --- |
+|  [L](./vector.choleskydecomposition.l.md) | <code>Matrix&lt;S&gt;</code> |  |
+
diff --git a/docs/api/vector.md b/docs/api/vector.md
index e556844..2e5269f 100644
--- a/docs/api/vector.md
+++ b/docs/api/vector.md
@@ -32,11 +32,13 @@
 |  [addRowToRow(matrix, targetRow, rowToAdd)](./vector.addrowtorow.md) | An elementary row operations which returns a new matrix whose row at <code>targetRow</code> has had the row at <code>rowToAdd</code> added to it. |
 |  [addScalarMultipleOfRowToRow(matrix, targetRow, rowToAdd, scalar)](./vector.addscalarmultipleofrowtorow.md) | An elementary row operations which returns a new matrix whose row at <code>targetRow</code> has had a scalar multiple of <code>rowToAdd</code> added to it. |
 |  [backwardDifferenceMatrix(binCount)](./vector.backwarddifferencematrix.md) | Builds a matrix that transforms a vector to a vector of backward differences |
+|  [calculateCholeskyDecomposition(A)](./vector.calculatecholeskydecomposition.md) | Uses the serial version of the Cholesky algorith to calculate the Cholesky decomposition of a matrix <code>A</code>. |
 |  [calculateEigenvalues(A, numIterations, throwOnFailure)](./vector.calculateeigenvalues.md) | Uses the QR algorithm to compute the eigenvalues of a matrix <code>A</code> |
 |  [calculateGeneralLeastSquares(dataPoints, functionTemplate, numberOfTerms)](./vector.calculategeneralleastsquares.md) | Calculates a regression model for an arbitrary function. |
 |  [calculateLinearLeastSquares(dataPoints)](./vector.calculatelinearleastsquares.md) | Calculates a linear regression model for the provided <code>dataPoints</code>. |
 |  [calculateLUDecomposition(A)](./vector.calculateludecomposition.md) | Uses the Doolittle algorithm to calculate the LU Decomposition of a matrix A. |
 |  [calculateQRDecomposition(A)](./vector.calculateqrdecomposition.md) | Uses the Graham-Schmidt process to calculate the QR decomposition of the matrix A. |
+|  [calculateSingularValueDecomposition(A)](./vector.calculatesingularvaluedecomposition.md) | Uses the Power Method to calculate the Singular Value Decomposition of a matrix <code>A</code> |
 |  [centralDifferenceMatrix(binCount)](./vector.centraldifferencematrix.md) | Builds a matrix that transforms a vector to a vector of central differences |
 |  [columnSumSupremumNorm(A)](./vector.columnsumsupremumnorm.md) | Calculates the 1-Norm of a matrix <code>A</code> |
 |  [crossProduct(first, second)](./vector.crossproduct.md) | Calculates the cross-product (vector-product) of two vectors. This is defined only for vectors with three dimensions. |
@@ -67,6 +69,7 @@
 |  [prettyPrint(num)](./vector.prettyprint.md) | Returns an easy-to-read string representing a <code>number</code> |
 |  [prettyPrint(vector)](./vector.prettyprint_1.md) | Returns an easy-to-read string representing the contents of a [Vector](./vector.vector.md) |
 |  [prettyPrint(matrix)](./vector.prettyprint_2.md) | Returns an easy-to-read string representing the contents of a [Matrix](./vector.matrix.md) |
+|  [rank(matrix)](./vector.rank.md) | Calculates the rank of a matrix |
 |  [reducedRowEchelonForm(matrix)](./vector.reducedrowechelonform.md) | Uses Gauss-Jordan elimination with pivoting to convert a matrix to Reduced Row-Echelon Form (RREF) |
 |  [rowEchelonForm(matrix)](./vector.rowechelonform.md) | Uses Gauss-Jordan elimination with pivoting to convert a matrix to Row-Echelon Form (REF) |
 |  [rowSumSupremumNorm(A)](./vector.rowsumsupremumnorm.md) | Calculates the Infinity-Norm of a matrix <code>A</code> |
@@ -80,6 +83,7 @@
 
 |  Interface | Description |
 |  --- | --- |
+|  [CholeskyDecomposition](./vector.choleskydecomposition.md) | The result of a Cholesky Decomposition |
 |  [EigenPair](./vector.eigenpair.md) | An eigenvector and its corresponding eigenvalue |
 |  [LeastSquaresApproximation](./vector.leastsquaresapproximation.md) | The result of a least squares approximation. |
 |  [LinearTransformation](./vector.lineartransformation.md) | An abstract linear transformation between vectors of type <code>V</code> and vectors of type <code>U</code>. |
@@ -87,6 +91,7 @@
 |  [Matrix](./vector.matrix.md) | A generalized Matrix - one of the core data types |
 |  [QRDecomposition](./vector.qrdecomposition.md) | The result of a QR decomposition. |
 |  [RowOperationResult](./vector.rowoperationresult.md) | The result of a row operation (<code>result</code>), and the matrix that we multiply by the original matrix to yield that result (<code>operator</code>) |
+|  [SingularValueDecomposition](./vector.singularvaluedecomposition.md) | The result of a Singular Value Decomposition |
 |  [Vector](./vector.vector.md) | A generalized Vector - one of the core data types |
 
 ## Type Aliases
diff --git a/docs/api/vector.rank.md b/docs/api/vector.rank.md
new file mode 100644
index 0000000..ad7d80d
--- /dev/null
+++ b/docs/api/vector.rank.md
@@ -0,0 +1,28 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [rank](./vector.rank.md)
+
+## rank() function
+
+Calculates the rank of a matrix
+
+<b>Signature:</b>
+
+```typescript
+export declare function rank<S>(matrix: Matrix<S>): number;
+```
+
+## Parameters
+
+|  Parameter | Type | Description |
+|  --- | --- | --- |
+|  matrix | <code>Matrix&lt;S&gt;</code> | the matrix for which to determine the rank |
+
+<b>Returns:</b>
+
+`number`
+
+## Remarks
+
+The rank of a matrix A is the dimension of the vector space spanned by the columns of A. Equivalently, it is the number of pivot entries in the row-echelon form of A, or the number of nonzero rows in the row echelon form of A.
+
diff --git a/docs/api/vector.singularvaluedecomposition.md b/docs/api/vector.singularvaluedecomposition.md
new file mode 100644
index 0000000..39d64d6
--- /dev/null
+++ b/docs/api/vector.singularvaluedecomposition.md
@@ -0,0 +1,22 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [SingularValueDecomposition](./vector.singularvaluedecomposition.md)
+
+## SingularValueDecomposition interface
+
+The result of a Singular Value Decomposition
+
+<b>Signature:</b>
+
+```typescript
+export interface SingularValueDecomposition<S> 
+```
+
+## Properties
+
+|  Property | Type | Description |
+|  --- | --- | --- |
+|  [Sigma](./vector.singularvaluedecomposition.sigma.md) | <code>Matrix&lt;S&gt;</code> |  |
+|  [U](./vector.singularvaluedecomposition.u.md) | <code>Matrix&lt;S&gt;</code> |  |
+|  [V](./vector.singularvaluedecomposition.v.md) | <code>Matrix&lt;S&gt;</code> |  |
+
diff --git a/docs/api/vector.singularvaluedecomposition.sigma.md b/docs/api/vector.singularvaluedecomposition.sigma.md
new file mode 100644
index 0000000..1998dd4
--- /dev/null
+++ b/docs/api/vector.singularvaluedecomposition.sigma.md
@@ -0,0 +1,11 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [SingularValueDecomposition](./vector.singularvaluedecomposition.md) &gt; [Sigma](./vector.singularvaluedecomposition.sigma.md)
+
+## SingularValueDecomposition.Sigma property
+
+<b>Signature:</b>
+
+```typescript
+Sigma: Matrix<S>;
+```
diff --git a/docs/api/vector.singularvaluedecomposition.u.md b/docs/api/vector.singularvaluedecomposition.u.md
new file mode 100644
index 0000000..1c77ced
--- /dev/null
+++ b/docs/api/vector.singularvaluedecomposition.u.md
@@ -0,0 +1,11 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [SingularValueDecomposition](./vector.singularvaluedecomposition.md) &gt; [U](./vector.singularvaluedecomposition.u.md)
+
+## SingularValueDecomposition.U property
+
+<b>Signature:</b>
+
+```typescript
+U: Matrix<S>;
+```
diff --git a/docs/api/vector.singularvaluedecomposition.v.md b/docs/api/vector.singularvaluedecomposition.v.md
new file mode 100644
index 0000000..35ccc06
--- /dev/null
+++ b/docs/api/vector.singularvaluedecomposition.v.md
@@ -0,0 +1,11 @@
+<!-- Do not edit this file. It is automatically generated by API Documenter. -->
+
+[Home](./index.md) &gt; [@josh-brown/vector](./vector.md) &gt; [SingularValueDecomposition](./vector.singularvaluedecomposition.md) &gt; [V](./vector.singularvaluedecomposition.v.md)
+
+## SingularValueDecomposition.V property
+
+<b>Signature:</b>
+
+```typescript
+V: Matrix<S>;
+```
diff --git a/package.json b/package.json
index 27d5869..3b4c60a 100644
--- a/package.json
+++ b/package.json
@@ -1,7 +1,7 @@
 {
   "name": "@josh-brown/vector",
   "author": "Joshua Brown <jbrown1.618@gmail.com>",
-  "version": "0.1.0",
+  "version": "0.2.0",
   "description": "A linear algebra library written in TypeScript that focuses on generality, extensibility, and ease of use.",
   "license": "ISC",
   "repository": {
diff --git a/src/algorithms/CholeskyDecomposition.ts b/src/algorithms/CholeskyDecomposition.ts
index 4cbc7fc..2e5775d 100644
--- a/src/algorithms/CholeskyDecomposition.ts
+++ b/src/algorithms/CholeskyDecomposition.ts
@@ -1,6 +1,10 @@
 import { Matrix } from '../types/matrix/Matrix';
 import { isHermitian } from '../utilities/MatrixProperties';
 
+/**
+ * The result of a Cholesky Decomposition
+ * @public
+ */
 export interface CholeskyDecomposition<S> {
   L: Matrix<S>;
 }
@@ -15,6 +19,7 @@ export interface CholeskyDecomposition<S> {
  *
  * A Cholesky decomposition only exists if `A` is symmetrix and positive-definite.
  * @param A - The matrix to decompose
+ * @public
  */
 export function calculateCholeskyDecomposition<S>(
   A: Matrix<S>
diff --git a/src/algorithms/GaussJordan.ts b/src/algorithms/GaussJordan.ts
index 0c012b1..4c85c7a 100644
--- a/src/algorithms/GaussJordan.ts
+++ b/src/algorithms/GaussJordan.ts
@@ -52,6 +52,24 @@ export function inverse<S>(matrix: Matrix<S>): Matrix<S> | undefined {
   }
 }
 
+/**
+ * Calculates the rank of a matrix
+ *
+ * @remarks
+ * The rank of a matrix A is the dimension of the vector space spanned by the columns of A.
+ * Equivalently, it is the number of pivot entries in the row-echelon form of A, or the number
+ * of nonzero rows in the row echelon form of A.
+ *
+ * @param matrix - the matrix for which to determine the rank
+ * @public
+ */
+export function rank<S>(matrix: Matrix<S>): number {
+  const zeroRow = matrix.vectorBuilder().zeros(matrix.getNumberOfColumns());
+  const ref = rowEchelonForm(matrix);
+  const nonZeroRows = ref.getRowVectors().filter(v => !v.equals(zeroRow));
+  return nonZeroRows.length;
+}
+
 /**
  * Uses Gauss-Jordan elimination with pivoting to convert a matrix to Reduced Row-Echelon Form (RREF)
  *
diff --git a/src/algorithms/SingularValueDecomposition.spec.ts b/src/algorithms/SingularValueDecomposition.spec.ts
new file mode 100644
index 0000000..667f69c
--- /dev/null
+++ b/src/algorithms/SingularValueDecomposition.spec.ts
@@ -0,0 +1,48 @@
+import { expect } from 'chai';
+import { describe, it } from 'mocha';
+import { NumberMatrix } from '../types/matrix/NumberMatrix';
+import { calculateSingularValueDecomposition } from './SingularValueDecomposition';
+import { NumberVector } from '../types/vector/NumberVector';
+
+describe('SingularValueDecomposition', () => {
+  const matrixBuilder = NumberMatrix.builder();
+  const vectorBuilder = NumberVector.builder();
+
+  describe('calculateSingularValueDecomposition', () => {
+    it('calculates the Singular Value Decomposition of a matrix A', () => {
+      const A = matrixBuilder.fromArray([[3, 2, 2], [2, 3, -2]]);
+
+      const oneOverRootTwo = 1 / Math.sqrt(2);
+      const oneOverRootEighteen = 1 / Math.sqrt(18);
+
+      const expectedSigma = matrixBuilder.diagonal(vectorBuilder.fromArray([5, 3]));
+      const expectedU = matrixBuilder.fromArray([
+        [oneOverRootTwo, oneOverRootTwo],
+        [oneOverRootTwo, -1 * oneOverRootTwo]
+      ]);
+      const expectedV = matrixBuilder.fromArray([
+        [oneOverRootTwo, oneOverRootEighteen],
+        [oneOverRootTwo, -1 * oneOverRootEighteen],
+        [0, 4 * oneOverRootEighteen]
+      ]);
+
+      const { U, Sigma, V } = calculateSingularValueDecomposition(A);
+
+      // Check that USigmaV* equals A
+      const USigmaVStar = U.multiply(Sigma.multiply(V.adjoint()));
+      expect(USigmaVStar.equals(A)).to.be.true;
+
+      expect(U.equals(expectedU)).to.be.true;
+      expect(V.equals(expectedV)).to.be.true;
+      expect(Sigma.equals(expectedSigma)).to.be.true;
+    });
+
+    it('calculates the SVD for a random 20x20 matrix', () => {
+      const A = matrixBuilder.random(20);
+      const { U, Sigma, V } = calculateSingularValueDecomposition(A);
+      // Check that USigmaV* equals A
+      const USigmaVStar = U.multiply(Sigma.multiply(V.adjoint()));
+      expect(USigmaVStar.equals(A)).to.be.true;
+    });
+  });
+});
diff --git a/src/algorithms/SingularValueDecomposition.ts b/src/algorithms/SingularValueDecomposition.ts
new file mode 100644
index 0000000..5b417c1
--- /dev/null
+++ b/src/algorithms/SingularValueDecomposition.ts
@@ -0,0 +1,90 @@
+import { Matrix } from '../types/matrix/Matrix';
+import { Vector } from '../types/vector/Vector';
+import { normalize, euclideanNorm } from './Norms';
+import { rank } from './GaussJordan';
+
+/**
+ * The result of a Singular Value Decomposition
+ * @public
+ */
+export interface SingularValueDecomposition<S> {
+  U: Matrix<S>;
+  Sigma: Matrix<S>;
+  V: Matrix<S>;
+}
+
+interface SingularValueAndVectors<S> {
+  u: Vector<S>;
+  sigma: S;
+  v: Vector<S>;
+}
+
+/**
+ * Uses the Power Method to calculate the Singular Value Decomposition of a matrix `A`
+ *
+ * @remarks
+ * A Singular Value Decomposition consists of orthogonal matrices `U` and `V`
+ * and a diagonal matrix `Sigma` such that _USigmaV* = A_
+ *
+ * @param A - the matrix to decompose
+ * @public
+ */
+export function calculateSingularValueDecomposition<S>(
+  A: Matrix<S>
+): SingularValueDecomposition<S> {
+  const ops = A.ops();
+  const builder = A.builder();
+  const vectorBuilder = A.vectorBuilder();
+  const r = rank(A);
+  const uColumns: Vector<S>[] = [];
+  const vColumns: Vector<S>[] = [];
+  const singularValues: S[] = [];
+
+  for (let i = 0; i < r; i++) {
+    const { u, sigma, v } = getFirstSingularValue(A);
+    uColumns.push(u);
+    singularValues.push(sigma);
+    vColumns.push(v);
+
+    const delta = u
+      .outerProduct(v)
+      .scalarMultiply(sigma)
+      .scalarMultiply(ops.negativeOne());
+    A = A.add(delta);
+  }
+
+  const U = builder.fromColumnVectors(uColumns);
+  const Sigma = builder.diagonal(vectorBuilder.fromArray(singularValues));
+  const V = builder.fromColumnVectors(vColumns);
+
+  return { U, Sigma, V };
+}
+
+function getFirstSingularValue<S>(A: Matrix<S>): SingularValueAndVectors<S> {
+  const ops = A.ops();
+  const vectorBuilder = A.vectorBuilder();
+
+  const x0 = vectorBuilder.ones(A.getNumberOfColumns());
+  const x = applyBToTheK(A, x0, 15); // TODO - what should k be?
+
+  const v = normalize(x);
+  if (v === undefined) throw new Error('TODO - Error');
+
+  const Av = A.apply(v);
+  const sigma = ops.fromNumber(euclideanNorm(Av));
+  const sigmaInverse = ops.getMultiplicativeInverse(sigma);
+  if (sigmaInverse === undefined) throw new Error('TODO - Error');
+  const u = Av.scalarMultiply(sigmaInverse);
+
+  return { u, sigma, v };
+}
+
+function applyBToTheK<S>(A: Matrix<S>, x0: Vector<S>, k: number): Vector<S> {
+  const At = A.adjoint();
+  let x = x0;
+  for (let i = 0; i < k; i++) {
+    x = A.apply(x);
+    x = At.apply(x);
+  }
+  return x;
+}
diff --git a/src/index.ts b/src/index.ts
index 83abbd1..dec90ed 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -1,3 +1,4 @@
+export * from './algorithms/CholeskyDecomposition';
 export * from './algorithms/CrossProduct';
 export * from './algorithms/Determinant';
 export * from './algorithms/Eigenvalues';
@@ -9,6 +10,7 @@ export * from './algorithms/LUDecomposition';
 export * from './algorithms/Norms';
 export * from './algorithms/QRDecomposition';
 export * from './algorithms/RowOperations';
+export * from './algorithms/SingularValueDecomposition';
 export * from './types/matrix/ArrayMatrix';
 export * from './types/matrix/ComplexMatrix';
 export * from './types/matrix/LinearTransformation';