# Optimal matrix chain multiplication in Java

Preliminaries

Given two matrices $$A = \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1q} \\ a_{21} & a_{22} & \dots & a_{2q} \\ \vdots & \vdots & \ddots & \vdots \\ a_{p1} & a_{p2} & \dots & a_{pq} \\ \end{pmatrix} \text{ and } B = \begin{pmatrix} b_{11} & b_{12} & \dots & b_{1r} \\ b_{21} & b_{22} & \dots & b_{2r} \\ \vdots & \vdots & \ddots & \vdots \\ b_{q1} & b_{q2} & \dots & b_{qr} \\ \end{pmatrix},$$ the product $M = AB$ is the matrix $$\begin{pmatrix} m_{11} & m_{12} & \dots & m_{1r} \\ m_{21} & m_{22} & \dots & m_{2r} \\ \vdots & \vdots & \ddots & \vdots \\ m_{p1} & m_{p2} & \dots & m_{pr} \\ \end{pmatrix},$$ where $$m_{ij} = \sum_{i = 1}^q a_{iq} b_{qj}.$$ The matrix product $AB$ is defined only if the number of columns in $A$ equals the number of rows in $B$. Also, multiplying the two matrices requires $pqr$ scalar multiplications.

Also, what comes to matrices, generally, two matrices $A$ and $B$ do not commute, i.e., $AB$ and $BA$ are not guaranteed to be the same even if both products are well-defined. However, matrices are associative: $A(BC) = (AB)C$.

Motivation

Suppose we are given a chain of three matrices $\langle A_1 \in \mathbb{R}^{10 \times 100}, A_2 \in \mathbb{R}^{100 \times 5}, A_3 \in \mathbb{R}^{5 \times 50} \rangle$. The parenthesation $A_1(A_2 A_3)$ produces $$(100 \times 5 \times 50) + (10 \times 100 \times 50) = 25000 + 50000 = 75000,$$ whereas the parenthesation $(A_1A_2)A_3$ produces $$(10 \times 100 \times 5) + (10 \times 5 \times 50) = 5000 + 2500 = 7500.$$ Clearly, the order of matrix multiplications in a matrix chain matters, and this post is about finding the optimal parenthesation that implies the least possible number of scalar multiplications and performing that optimal sequence of matrix multiplications.

MatrixChainMultiplier.java:

package net.coderodde.matrix;

/**
* This abstract class defines the API for matrix chain multiplication
* algorithms.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Dec 19, 2015)
*/
public interface MatrixChainMultiplier {

/**
* Computes the product of the matrix chain.
*
* @param matrices an array of matrices to multiply.
* @return         the matrix chain product.
*/
public Matrix multiply(Matrix... matrices);
}

OptimalMatrixChainMultiplier.java:

package net.coderodde.matrix.support;

import net.coderodde.matrix.Matrix;
import net.coderodde.matrix.MatrixChainMultiplier;

/**
* This class implements an optimal matrix chain multiplying algorithm, which
* performs the minimum number of scalar multiplications. This algorithm is
* described in the 3rd edition of "Introduction to Algorithms", Chapter 15:
* Dynamic Programming, page 375.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Dec 19, 2015)
*/
public class OptimalMatrixChainMultiplier implements MatrixChainMultiplier {

@Override
public Matrix multiply(Matrix... matrices) {
switch (matrices.length) {
case 0: {
return null;
}

case 1: {
return matrices[0].clone();
}
}

Matrix[] cloneMatrices = new Matrix[matrices.length];

for (int i = 0; i < matrices.length; ++i) {
cloneMatrices[i] = matrices[i].clone();
}

return computeOptimalParenthesation(cloneMatrices);
}

private static int[] extractDimensionArray(Matrix[] chain) {
int[] input = new int[chain.length + 1];

for (int i = 0; i < chain.length; ++i) {
input[i] = chain[i].getHeight();
}

input[input.length - 1] = chain[chain.length - 1].getWidth();
return input;
}

private Matrix computeOptimalParenthesation(Matrix[] matrices) {
int[] p = extractDimensionArray(matrices);
int n = p.length - 1;
int[][] m = new int[n + 1][n + 1];
int[][] s = new int[n + 1][n + 1];

for (int len = 2; len <= n; ++len) {
for (int i = 1; i <= n - len + 1; ++i) {
int j = i + len - 1;
m[i][j] = Integer.MAX_VALUE;

for (int k = i; k <= j - 1; ++k) {
int q = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];

if (m[i][j] > q) {
m[i][j] = q;
s[i][j] = k;
}
}
}
}

return process(s, matrices, 1, n);
}

private Matrix process(int[][] s, Matrix[] matrices, int i, int j) {
if (i == j) {
return matrices[i - 1];
}

return process(s, matrices, i, s[i][j]).multiply(process(s,
matrices,
s[i][j] + 1,
j));
}
}

NaiveMatrixChainMultiplier.java:

package net.coderodde.matrix.support;

import net.coderodde.matrix.IncompatibleMatrixException;
import net.coderodde.matrix.Matrix;
import net.coderodde.matrix.MatrixChainMultiplier;

/**
* This class implements a naive algorithm for multiplying matrix chains.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Dec 19, 2015)
*/
public class NaiveMatrixChainMultiplier implements MatrixChainMultiplier {

@Override
public Matrix multiply(Matrix... matrices) {
switch (matrices.length) {
case 0: {
return null;
}

case 1: {
return matrices[0].clone();
}
}

Matrix current = matrices[0].clone();

for (int i = 1; i < matrices.length; ++i) {
current = current.multiply(matrices[i]);
}

return current;
}
}

Matrix.java:

package net.coderodde.matrix;

import java.util.Arrays;

/**
* This class implements a matrix data type.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Dec 19, 2015)
*/
public class Matrix implements Cloneable {

private final double[][] matrix;

/**
* Creates a square matrix with width and height equal to {@code dimension}.
*
* @param dimension the width and height of the square matrix.
*/
public Matrix(int dimension) {
this(dimension, dimension);
}

public Matrix(int width, int height) {
checkWidth(width);
checkHeight(height);
this.matrix = new double[height][width];
}

public int getWidth() {
return matrix[0].length;
}

public int getHeight() {
return matrix.length;
}

public double read(int x, int y) {
return matrix[y][x];
}

public void write(int x, int y, double value) {
matrix[y][x] = value;
}

public Matrix multiply(Matrix rightHandMatrix) {
checkDimensions(rightHandMatrix);

int thisWidth = getWidth();
int thisHeight = getHeight();
int otherWidth = rightHandMatrix.getWidth();

Matrix ret = new Matrix(otherWidth, thisHeight);

for (int thisY = 0; thisY < thisHeight; ++thisY) {
for (int otherX = 0; otherX < otherWidth; ++otherX) {
double sum = 0.0;

for (int i = 0; i < thisWidth; ++i) {
}

ret.write(otherX, thisY, sum);
}
}

return ret;
}

public Matrix clone() {
Matrix ret = new Matrix(getWidth(), getHeight());
int width = getWidth();
int height = getHeight();

for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
}
}

return ret;
}

@Override
public boolean equals(Object o) {
if (o == null) {
return false;
}

if (!getClass().equals(o.getClass())) {
return false;
}

Matrix other = (Matrix) o;

return Arrays.deepEquals(matrix, other.matrix);
}

@Override
public int hashCode() {
return Arrays.deepHashCode(matrix);
}

@Override
public String toString() {
StringBuilder sb = new StringBuilder();

int width = getWidth();
int height = getHeight();
int maxIntegerPartLength = 0;
int maxDecimalPartLength = 0;

for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
String[] parts = s.split("\\.");

if (maxIntegerPartLength < parts[0].length()) {
maxIntegerPartLength = parts[0].length();
}

if (maxDecimalPartLength < parts[1].length()) {
maxDecimalPartLength = parts[1].length();
}
}
}

String fmt = "%" + (maxIntegerPartLength + maxDecimalPartLength + 1)
+ "." + maxDecimalPartLength + "f";

for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {

if (x < width - 1) {
sb.append(' ');
}
}

if (y < height - 1) {
sb.append('\n');
}
}

return sb.toString();
}

private void checkWidth(int width) {
checkIsPositive(width);
}

private void checkHeight(int height) {
checkIsPositive(height);
}

private void checkDimensions(Matrix rightHandMatrix) {
if (this.getWidth() != rightHandMatrix.getHeight()) {
throw new IncompatibleMatrixException(
"Dimension mismatch. The number of columns in this " +
"matrix (" + getWidth() + ") does not equal the number " +
"of rows in the right hand matrix (" +
rightHandMatrix.getHeight() + ").");
}
}

private void checkIsPositive(int number) {
if (number < 1) {
throw new IllegalArgumentException(
"The input number is not positive. Received: " + number);
}
}
}

IncompatibleMatrixException.java:

package net.coderodde.matrix;

/**
* This class implements an exception that is thrown whenever multiplying two
* <b>incompatible</b> matrices. Whenever computing the product {@code AB}, the
* two matrices {@code A} and {@code B} are incompatible if and only if the
* number of columns in {@code A} does not equal the number of rows of
* {@code B}.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Dec 19, 2015)
*/
public class IncompatibleMatrixException extends RuntimeException {

public IncompatibleMatrixException(String message) {
super(message);
}
}

PerformanceDemo.java:

package net.coderodde.matrix;

import java.util.Random;
import net.coderodde.matrix.support.NaiveMatrixChainMultiplier;
import net.coderodde.matrix.support.OptimalMatrixChainMultiplier;

public class PerformanceDemo {

private static final int MATRIX_CHAIN_LENGTH = 6;
private static final int MAXIMUM_WIDTH = 500;
private static final int MAXIMUM_HEIGHT = 500;
private static final int MINIMUM_VALUE = -5;
private static final int MAXIMUM_VALUE = 5;
private static final double E = 0.001;

public static void main(final String... args) {
long seed = System.nanoTime();
Random random = new Random(seed);
Matrix[] chain = createRandomMatrixChain(MATRIX_CHAIN_LENGTH,
MAXIMUM_WIDTH,
MAXIMUM_HEIGHT,
MINIMUM_VALUE,
MAXIMUM_VALUE,
random);
System.out.println("Seed = " + seed);

Matrix result1 = profile(new NaiveMatrixChainMultiplier(), chain);
Matrix result2 = profile(new OptimalMatrixChainMultiplier(), chain);

System.out.println("Matrices equal: " + equals(result1, result2, E));
}

private static boolean equals(Matrix m1, Matrix m2, double e) {
if (m1.getWidth() != m2.getWidth()) {
return false;
}

if (m1.getHeight() != m2.getHeight()) {
return false;
}

for (int y = 0; y < m1.getHeight(); ++y) {
for (int x = 0; x < m1.getWidth(); ++x) {
return false;
}
}
}

return true;
}

private static Matrix profile(MatrixChainMultiplier multiplier,
Matrix[] chain) {
long startTime = System.nanoTime();
Matrix ret = multiplier.multiply(chain);
long endTime = System.nanoTime();

System.out.println(multiplier.getClass().getSimpleName() + " in " +
(int)((endTime - startTime) / 1e6) +
" milliseconds.");

return ret;
}

private static Matrix[] createRandomMatrixChain(int matrices,
int maximumWidth,
int maximumHeight,
int minimumValue,
int maximumValue,
Random random) {
Matrix[] ret = new Matrix[matrices];

if (ret.length == 0) {
return ret;
}

ret[0] = createRandomMatrix(randomInt(1, maximumWidth, random),
randomInt(1, maximumHeight, random),
minimumValue,
maximumValue,
random);

int previousWidth = ret[0].getWidth();

for (int i = 1; i < matrices; ++i) {
ret[i] = createRandomMatrix(randomInt(1, maximumWidth, random),
previousWidth,
minimumValue,
maximumValue,
random);
previousWidth = ret[i].getWidth();
}

return  ret;
}

private static Matrix createRandomMatrix(int width,
int height,
int minimumValue,
int maximumValue,
Random random) {
Matrix ret = new Matrix(width, height);

for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
ret.write(x, y, randomInt(minimumValue, maximumValue, random));
}
}

return ret;
}

private static int randomInt(int minimumValue,
int maximumValue,
Random random) {
return random.nextInt(maximumValue - minimumValue + 1) + minimumValue;
}
}

MatrixTest.java:

package net.coderodde.matrix;

import org.junit.Test;
import static org.junit.Assert.*;

public class MatrixTest {

private static final double E = 1e-4;
private Matrix m;

@Test
public void testSmallestMatrixDoesNotThrow() {
new Matrix(1, 1);
}

@Test(expected = IllegalArgumentException.class)
public void testMatrixThrowsOnNonPositiveWidth() {
new Matrix(0, 1);
}

@Test(expected = IllegalArgumentException.class)
public void testMatrixThrowsOnNonPositiveHeight() {
new Matrix(1, 0);
}

@Test(expected = IllegalArgumentException.class)
public void testMatrixThrowsOnNonPositiveWidthAndHeight() {
new Matrix(0, 0);
}

@Test(expected = IncompatibleMatrixException.class)
public void testMatrixThrowsOnMultiplyintIncompatibleMatrix() {
new Matrix(3, 2).multiply(new Matrix(2, 2));
}

@Test
public void testGetWidth() {
m = new Matrix(2, 3);
assertEquals(2, m.getWidth());

m = new Matrix(5, 3);
assertEquals(5, m.getWidth());

m = new Matrix(5, 4);
assertEquals(5, m.getWidth());
}

@Test
public void testGetHeight() {
m = new Matrix(4, 6);
assertEquals(6, m.getHeight());

m = new Matrix(2, 6);
assertEquals(6, m.getHeight());

m = new Matrix(2, 5);
assertEquals(5, m.getHeight());
}

@Test
m = new Matrix(3, 3);

for (int y = 0; y < m.getHeight(); ++y) {
for (int x = 0; x < m.getWidth(); ++x) {
}
}

m.write(1, 1, 4.0);

m.write(2, 2, 5.0);

m.write(2, 2, 9.5);

m.write(1, 1, 2.0);
}

@Test
public void testWrite() {
m = new Matrix(4, 3);

m.write(3, 1, 4.6);

m.write(3, 1, 2.6);
}

@Test
public void testMultiply() {
// 1, 2,
// 3, 4
// 5, 6
m = new Matrix(2, 3);

int num = 1;

for (int y = 0; y < m.getHeight(); ++y) {
for (int x = 0; x < m.getWidth(); ++x) {
m.write(x, y, num++);
}
}

// 3,  2
// 1, -1

Matrix m2 = new Matrix(2, 2);

m2.write(0, 0, 3);
m2.write(1, 0, 2);
m2.write(0, 1, 1);
m2.write(1, 1, -1);

// 3, 2
// 1, -1
// 1, 2              5,  0
// 3, 4     x      = 13, 2
// 5, 6              21, 4
Matrix result = m.multiply(m2);

assertEquals(2, result.getWidth());
assertEquals(3, result.getHeight());

}

@Test
public void testEquals() {
m = new Matrix(2, 2);
Matrix m2 = new Matrix(2, 3);

assertFalse(m.equals(m2));

m2 = new Matrix(3, 2);
assertFalse(m.equals(m2));

m2 = new Matrix(2, 2);
assertEquals(m, m2);

m.write(1, 0, 3.14);
assertFalse(m.equals(m2));

m2.write(1, 0, 3.14);
assertEquals(m, m2);

m2.write(0, 1, 2.78);
assertFalse(m.equals(m2));

assertEquals(m, m);
assertEquals(m2, m2);

assertFalse(m.equals(null));
assertFalse(m.equals(new Object()));
}

@Test
public void testHashCode() {
m = new Matrix(3, 3);
m.write(1, 1, 5.0);

Matrix m2 = new Matrix(3, 3);

assertFalse(m.hashCode() == m2.hashCode());

m2.write(1, 1, 4.0);
assertFalse(m.hashCode() == m2.hashCode());

m2.write(1, 1, 5.0);
assertEquals(m.hashCode(), m2.hashCode());
}

@Test
public void testToString() {
System.out.println("toString():");

m = new Matrix(3, 4);
m.write(2, 2, 3.14);
m.write(2, 3, 0.0002);
m.write(0, 1, 268.9);

System.out.println(m);
}
}

NaiveMatrixChainMultiplierTest.java:

package net.coderodde.matrix.support;

import net.coderodde.matrix.Matrix;
import org.junit.Test;
import static org.junit.Assert.*;

public class NaiveMatrixChainMultiplierTest {

private static final double E = 1e-5;

@Test
public void testMultiply() {
Matrix[] matrices = new Matrix[3];

matrices[0] = new Matrix(2, 2);
matrices[1] = new Matrix(3, 2);
matrices[2] = new Matrix(2, 3);

int num = 1;

for (Matrix m : matrices) {
int width = m.getWidth();
int height = m.getHeight();

for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
m.write(x, y, num++);
}
}
}

Matrix result = new NaiveMatrixChainMultiplier().multiply(matrices);

assertEquals(2, result.getWidth());
assertEquals(2, result.getHeight());

}
}

OptimalMatrixChainMultiplier.java:

package net.coderodde.matrix.support;

import net.coderodde.matrix.Matrix;
import org.junit.Test;
import static org.junit.Assert.*;

public class OptimalMatrixChainMultiplierTest {

private static final double E = 1e-5;

@Test
public void testMultiply() {
Matrix[] matrices = new Matrix[3];

matrices[0] = new Matrix(2, 2);
matrices[1] = new Matrix(3, 2);
matrices[2] = new Matrix(2, 3);

int num = 1;

for (Matrix m : matrices) {
int width = m.getWidth();
int height = m.getHeight();

for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
m.write(x, y, num++);
}
}
}

Matrix result = new OptimalMatrixChainMultiplier().multiply(matrices);

assertEquals(2, result.getWidth());
assertEquals(2, result.getHeight());

}
}

Performance figures

Seed = 277908980432723
NaiveMatrixChainMultiplier in 420 milliseconds.
OptimalMatrixChainMultiplier in 27 milliseconds.
Matrices equal: true

Note

If you increase the constants provided in the demonstration, you may get the message that the two result matrices are not equal. This happens due to the fact that entries in the result matrices may get really large and a slightest difference in the mantissas of two corresponding entries will imply large absolute difference.

Critique request

Is my naming conventions/API design/algorithm implementation in order? Anything else to improve?

Your tests are very repetitive, largely copy-pasted may I dare to say. Copy-pasting is not an acceptable coding tool.

I suggest you merge your OptimalMatrixChainMultiplierTest.java and NaiveMatrixChainMultiplierTest.java: in a single MatrixChainMultiplierTest.java organized this way:

import ...

public class MatrixChainMultiplierTest {
private static Matrix InputMatrix() {
// Code to generate the matrix goes here.
}

@Test
public void testMultiplyNaive() {
Matrix result = new OptimalMatrixChainMultiplier().multiply(InputMatrix());

assertEquals(2, result.getWidth());
assertEquals(2, result.getHeight());

}

@Test
public void testMultiplyNaive() {
Matrix result = new NaiveMatrixChainMultiplier().multiply(InputMatrix());

assertEquals(2, result.getWidth());
assertEquals(2, result.getHeight());

}

There is still some repetition, but this already has 15 less duplicate lines than your version.

• The first test should be called testMultiplyOptimal. Oct 31 '19 at 16:29

(not intending to delve into demo&test.)

naming conventions/API design

(Are you aware of JAMA?)

• The interface consist of a single function, documented and appropriately named in the usual mix of cases - formally exemplary, but don't call it abstract class. Food for thought (think package summary):

/** Performs chain multiplication of {@code Matrix} to allow
*   improvements on {@code M1.multiply(M2)...multiply(Mn)}. */

(The following may not really be part of your question - ?)

• Matrix:

• Consider defining an interface. I'd go with get(row, col) and set(row, col, val). (I presume other operations have been left out to avoid clutter irrelevant to chained multiplication.)
• For Matrix.clone() to be a remedy against modification using the return value, it should state

/** Creates and returns a copy of this Matrix.
*  @return a {@code Matrix} not identical, but equal to
*   {@code this} with no means to modify {@code this} */

(BTW, not using super.clone() isn't illegal, but easily breaks expectations. See documentation of Object.clone() for details.)

• equals()' use of getClass().equals() is more restrictive than

return o instanceof Matrix
&& Arrays.deepEquals(matrix, ((Matrix)o).matrix);

• IncompatibleMatrixException: not providing a default constructor makes not doing the right thing (providing an illuminating message) look silly - nice. Consider extending IllegalArgumentException.

algorithm implementation

While documentation of the interface has priority, I'm not happy with a text book reference as the only comment.
In both MatrixChainMultiplier implementations, I don't see the need to clone the matrices. I think this is a documentation issue: 1) Matrix.multiply(right) doesn't specify anything, in particular, neither whether it modifies this or right nor whether it might return either one; 2) MatrixChainMultiplier.multiply() "extends Matrix.multiply".
Aware of the presentation of the algorithm in CLRS, I'm still against single letter variable names (with qualified exceptions for loop control variables of innermost loops (and n for the count nobody could mistake another for)). Is that s for splits and m for minimums (as far as evaluated)? I don't like guessing. (Speaking of which: the choice of int for the minimums deserves a comment. I'd use
long ops = (long)min[i][k] + min[k+1][j] + dim[i-1] * (long)dim[k] * dim[j];
(long q = (long)m[i][k] + m[k + 1][j] + p[i - 1] * (long)p[k] * p[j];) - just in case.)

(The uncomely business of parenthesization being less uncommon than parenthesation with parenthesisation in the middle makes me suggest association. Reminds of the origin of this degree of freedom, and optimalGrouping is dull.) computeOptimalParenthesation(matrices) should do just that, returning the two-dimensional array of splits (to the invoking multiply() to invoke process() with). Refactoring

private int[][] computeOptimalParenthesation(Matrix[] matrices) {
int[] dimensions = extractDimensionArray(matrices);
return computeOptimalParenthesation(dimensions);
}

highlights an opportunity for result reuse: if the sequence of dimensions has been evaluated in a past invocation, the optimal association should be the same this time.

Zhu, Wu, & Wang claim in a 2014 JofCIS article that their "Efficient Algorithm for Chain Multiplication of Dense Matrices" is both simpler and faster than Hu/Shins 1981 O(nlogn) algorithm (is that real Pascal code?) - any takers for a free implementation?

• (Added a case 2: to OptimalMatrixChainMultiplier.multiply(): meh.) Jan 31 '16 at 9:35