Matrix chain multiplication

From Rosetta Code
Task
Matrix chain multiplication
You are encouraged to solve this task according to the task description, using any language you may know.
Problem

Using the most straightforward algorithm (which we assume here), computing the product of two matrices of dimensions (n1,n2) and (n2,n3) requires n1*n2*n3 FMA operations. The number of operations required to compute the product of matrices A1, A2... An depends on the order of matrix multiplications, hence on where parens are put. Remember that the matrix product is associative, but not commutative, hence only the parens can be moved.

For instance, with four matrices, one can compute A(B(CD)), A((BC)D), (AB)(CD), (A(BC))D, (AB)C)D. The number of different ways to put the parens is a Catalan number, and grows exponentially with the number of factors.

Here is an example of computation of the total cost, for matrices A(5,6), B(6,3), C(3,1):

  • AB costs 5*6*3=90 and produces a matrix of dimensions (5,3), then (AB)C costs 5*3*1=15. The total cost is 105.
  • BC costs 6*3*1=18 and produces a matrix of dimensions (6,1), then A(BC) costs 5*6*1=30. The total cost is 48.

In this case, computing (AB)C requires more than twice as many operations as A(BC). The difference can be much more dramatic in real cases.

Task

Write a function which, given a list of the successive dimensions of matrices A1, A2... An, of arbitrary length, returns the optimal way to compute the matrix product, and the total cost. Any sensible way to describe the optimal solution is accepted. The input list does not duplicate shared dimensions: for the previous example of matrices A,B,C, one will only pass the list [5,6,3,1] (and not [5,6,6,3,3,1]) to mean the matrix dimensions are respectively (5,6), (6,3) and (3,1). Hence, a product of n matrices is represented by a list of n+1 dimensions.

Try this function on the following two lists:

  • [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
  • [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]

To solve the task, it's possible, but not required, to write a function that enumerates all possible ways to parenthesize the product. This is not optimal because of the many duplicated computations, and this task is a classic application of dynamic programming.

See also Matrix chain multiplication on Wikipedia.

11l

Translation of: Nim
T Optimizer
   [Int] dims
   [[Int]] m, s

   F (dims)
      .dims = dims

   F findMatrixChainOrder()
      V n = .dims.len - 1
      .m = [[0] * n] * n
      .s = [[0] * n] * n

      L(lg) 1 .< n
         L(i) 0 .< n - lg
            V j = i + lg
            .m[i][j] = 7FFF'FFFF
            L(k) i .< j
               V cost = .m[i][k] + .m[k + 1][j] + .dims[i] * .dims[k + 1] * .dims[j + 1]
               I cost < .m[i][j]
                  .m[i][j] = cost
                  .s[i][j] = k

   F optimalChainOrder(i, j)
      I i == j
         R String(Char(code' i + ‘A’.code))
      E
         R ‘(’(.optimalChainOrder(i, .s[i][j]))‘’
            ‘’(.optimalChainOrder(.s[i][j] + 1, j))‘)’

V Dims1 = [5, 6, 3, 1]
V Dims2 = [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
V Dims3 = [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]

L(dims) [Dims1, Dims2, Dims3]
   V opt = Optimizer(dims)
   opt.findMatrixChainOrder()
   print(‘Dims:  ’dims)
   print(‘Order: ’opt.optimalChainOrder(0, dims.len - 2))
   print(‘Cost:  ’opt.m[0][dims.len - 2])
   print(‘’)
Output:
Dims:  [5, 6, 3, 1]
Order: (A(BC))
Cost:  48

Dims:  [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order: ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost:  38120

Dims:  [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order: (A((((((BC)D)(((EF)G)H))I)J)K))
Cost:  1773740

Ada

This example implements the pseudocode in the reference Wiki page. The pseudocode states that the index values for the array to multiply begin at 0 while the cost and order matrices employ index values beginning at 1. Ada supports this pseudocode directly because Ada allows the programmer to define the index range for any array type. This Ada example is implemented using a simple package and a main procedure. The package specification is:

package mat_chain is
   type Vector is array (Natural range <>) of Integer;
   procedure Chain_Multiplication (Dims : Vector);
end mat_chain;

The implementation or body of the package is:

with Ada.Text_IO;           use Ada.Text_IO;
with Ada.Strings.Unbounded; use Ada.Strings.Unbounded;

package body mat_chain is
   type Result_Matrix is
     array (Positive range <>, Positive range <>) of Integer;

   --------------------------
   -- Chain_Multiplication --
   --------------------------

   procedure Chain_Multiplication (Dims : Vector) is
      n : Natural := Dims'Length - 1;
      S : Result_Matrix (1 .. n, 1 .. n);
      m : Result_Matrix (1 .. n, 1 .. n);
      procedure Print (Item : Vector) is
      begin
         Put ("Array Dimension  = (");
         for I in Item'Range loop
            Put (Item (I)'Image);
            if I < Item'Last then
               Put (",");
            else
               Put (")");
            end if;
         end loop;
         New_Line;
      end Print;

      procedure Chain_Order (Item : Vector) is
         J    : Natural;
         Cost : Natural;
         Temp : Natural;

      begin
         for idx in 1 .. n loop
            m (idx, idx) := 0;
         end loop;

         for Len in 2 .. n loop
            for I in 1 .. n - Len + 1 loop
               J        := I + Len - 1;
               m (I, J) := Integer'Last;
               for K in I .. J - 1 loop
                  Temp := Item (I - 1) * Item (K) * Item (J);
                  Cost := m (I, K) + m (K + 1, J) + Temp;
                  if Cost < m (I, J) then
                     m (I, J) := Cost;
                     S (I, J) := K;
                  end if;
               end loop;
            end loop;
         end loop;
      end Chain_Order;

      function Optimal_Parens return String is
         function Construct
           (S : Result_Matrix; I : Natural; J : Natural)
            return Unbounded_String
         is
            Us         : Unbounded_String := Null_Unbounded_String;
            Char_Order : Character;
         begin
            if I = J then
               Char_Order := Character'Val (I + 64);
               Append (Source => Us, New_Item => Char_Order);
               return Us;
            else
               Append (Source => Us, New_Item => '(');
               Append (Source => Us, New_Item => Construct (S, I, S (I, J)));
               Append (Source => Us, New_Item => '*');
               Append
                 (Source => Us, New_Item => Construct (S, S (I, J) + 1, J));
               Append (Source => Us, New_Item => ')');
               return Us;
            end if;
         end Construct;

      begin
         return To_String (Construct (S, 1, n));

      end Optimal_Parens;

   begin
      Chain_Order (Dims);
      Print (Dims);
      Put_Line ("Cost             = " & Integer'Image (m (1, n)));
      Put_Line ("Optimal Multiply = " & Optimal_Parens);
   end Chain_Multiplication;

end mat_chain;

The main procedure is:

with Mat_Chain; use Mat_Chain;
with Ada.Text_IO; use Ada.Text_IO;

procedure chain_main is
   V1 : Vector := (5, 6, 3, 1);
   V2 : Vector := (1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2);
   V3 : Vector := (1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10);
begin
   Chain_Multiplication(V1);
   New_Line;
   Chain_Multiplication(V2);
   New_Line;
   Chain_Multiplication(V3);
end chain_main;
Output:
Array Dimension  = ( 5, 6, 3, 1)
Cost             =  48
Optimal Multiply = (A*(B*C))

Array Dimension  = ( 1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2)
Cost             =  38120
Optimal Multiply = ((((((((A*B)*C)*D)*E)*F)*G)*(H*(I*J)))*(K*L))

Array Dimension  = ( 1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10)
Cost             =  1773740
Optimal Multiply = (A*((((((B*C)*D)*(((E*F)*G)*H))*I)*J)*K))

C

Translation of: Kotlin
#include <stdio.h>
#include <limits.h>
#include <stdlib.h>

int **m;
int **s;

void optimal_matrix_chain_order(int *dims, int n) {
    int len, i, j, k, temp, cost;
    n--;
    m = (int **)malloc(n * sizeof(int *));
    for (i = 0; i < n; ++i) {
        m[i] = (int *)calloc(n, sizeof(int));
    }

    s = (int **)malloc(n * sizeof(int *));
    for (i = 0; i < n; ++i) {
        s[i] = (int *)calloc(n, sizeof(int));
    }

    for (len = 1; len < n; ++len) {
        for (i = 0; i < n - len; ++i) {
            j = i + len;
            m[i][j] = INT_MAX;
            for (k = i; k < j; ++k) {
                temp = dims[i] * dims[k + 1] * dims[j + 1];
                cost = m[i][k] + m[k + 1][j] + temp;
                if (cost < m[i][j]) {
                    m[i][j] = cost;
                    s[i][j] = k;
                }
            }
        }
    }
}

void print_optimal_chain_order(int i, int j) {
    if (i == j)
        printf("%c", i + 65);
    else {
        printf("(");
        print_optimal_chain_order(i, s[i][j]);
        print_optimal_chain_order(s[i][j] + 1, j);
        printf(")");
    }
}

int main() {
    int i, j, n;
    int a1[4]  = {5, 6, 3, 1};
    int a2[13] = {1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2};
    int a3[12] = {1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10};
    int *dims_list[3] = {a1, a2, a3};
    int sizes[3] = {4, 13, 12};
    for (i = 0; i < 3; ++i) {
        printf("Dims  : [");
        n = sizes[i];
        for (j = 0; j < n; ++j) {
            printf("%d", dims_list[i][j]);
            if (j < n - 1) printf(", "); else printf("]\n");
        }        
        optimal_matrix_chain_order(dims_list[i], n);
        printf("Order : ");
        print_optimal_chain_order(0, n - 2);
        printf("\nCost  : %d\n\n", m[0][n - 2]);
        for (j = 0; j <= n - 2; ++j) free(m[j]);
        free(m);
        for (j = 0; j <= n - 2; ++j) free(s[j]);
        free(s);
    }
    return 0;
}
Output:
Dims  : [5, 6, 3, 1]
Order : (A(BC))
Cost  : 48

Dims  : [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order : ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost  : 38120

Dims  : [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order : (A((((((BC)D)(((EF)G)H))I)J)K))
Cost  : 1773740

C#

Translation of: Kotlin
using System;

class MatrixChainOrderOptimizer {
    private int[,] m;
    private int[,] s;

    void OptimalMatrixChainOrder(int[] dims) {
        int n = dims.Length - 1;
        m = new int[n, n];
        s = new int[n, n];
        for (int len = 1; len < n; ++len) {
            for (int i = 0; i < n - len; ++i) {
                int j = i + len;
                m[i, j] = Int32.MaxValue;
                for (int k = i; k < j; ++k) {
                    int temp = dims[i] * dims[k + 1] * dims[j + 1];
                    int cost = m[i, k] + m[k + 1, j] + temp;
                    if (cost < m[i, j]) {
                        m[i, j] = cost;
                        s[i, j] = k;
                    }
                }
            }
        }
    }

    void PrintOptimalChainOrder(int i, int j) {
        if (i == j)
            Console.Write((char)(i + 65));
        else {
            Console.Write("(");
            PrintOptimalChainOrder(i, s[i, j]);
            PrintOptimalChainOrder(s[i, j] + 1, j);
            Console.Write(")");
        }
    }

    static void Main() {
        var mcoo = new MatrixChainOrderOptimizer();
        var dimsList = new int[3][];
        dimsList[0] = new int[4] {5, 6, 3, 1};
        dimsList[1] = new int[13] {1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2};
        dimsList[2] = new int[12] {1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10};
        for (int i = 0; i < dimsList.Length; ++i) {
            Console.Write("Dims  : [");
            int n = dimsList[i].Length;
            for (int j = 0; j < n; ++j) {
                Console.Write(dimsList[i][j]);
                if (j < n - 1)
                    Console.Write(", ");
                else
                    Console.WriteLine("]");
            }
            mcoo.OptimalMatrixChainOrder(dimsList[i]);
            Console.Write("Order : ");
            mcoo.PrintOptimalChainOrder(0, n - 2);
            Console.WriteLine("\nCost  : {0}\n",  mcoo.m[0, n - 2]);
        }
    }
}
Output:
Dims  : [5, 6, 3, 1]
Order : (A(BC))
Cost  : 48

Dims  : [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order : ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost  : 38120

Dims  : [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order : (A((((((BC)D)(((EF)G)H))I)J)K))
Cost  : 1773740

C++

#include <cstdint>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

constexpr int32_t MAXIMUM_VALUE = 2'147'483'647;

std::vector<std::vector<int32_t>> cost;
std::vector<std::vector<int32_t>> order;

void print_vector(const std::vector<int32_t>& list) {
	std::cout << "[";
	for ( uint64_t i = 0; i < list.size() - 1; ++i ) {
		std::cout << list[i] << ", ";
	}
	std::cout << list.back() << "]" << std::endl;
}

int32_t matrix_chain_order(const std::vector<int32_t>& dimensions) {
	const uint64_t size = dimensions.size() - 1;
	cost = { size, std::vector<int32_t>(size, 0) };
	order = { size, std::vector<int32_t>(size, 0) };

	for ( uint64_t m = 1; m < size; ++m ) {
		for ( uint64_t i = 0; i < size - m; ++i ) {
			int32_t j = i + m;
			cost[i][j] = MAXIMUM_VALUE;
			for ( int32_t k = i; k < j; ++k ) {
				int32_t current_cost = cost[i][k] + cost[k + 1][j]
					+ dimensions[i] * dimensions[k + 1] * dimensions[j + 1];
				if ( current_cost < cost[i][j] ) {
					cost[i][j] = current_cost;
					order[i][j] = k;
				}
			}
		}
	}
	return cost[0][size - 1];
}

std::string get_optimal_parenthesizations(const std::vector<std::vector<int32_t>>& order,
										  const uint64_t& i, const uint64_t& j) {
	if ( i == j ) {
		std::string result(1, char(i + 65));
		return result;
	} else {
		std::stringstream stream;
		stream << "(" << get_optimal_parenthesizations(order, i, order[i][j])
			   << " * " << get_optimal_parenthesizations(order, order[i][j] + 1, j) << ")";
		return stream.str();
	}
}

void matrix_chain_multiplication(const std::vector<int32_t>& dimensions) {
	std::cout << "Array Dimension  = "; print_vector(dimensions);
	std::cout << "Cost             = " << matrix_chain_order(dimensions) << std::endl;
	std::cout << "Optimal Multiply = "
			  <<  get_optimal_parenthesizations(order, 0, order.size() - 1) << std::endl << std::endl;
}

int main() {
	matrix_chain_multiplication({ 5, 6, 3, 1 });
	matrix_chain_multiplication({ 1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2 });
	matrix_chain_multiplication({ 1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10 });
}
Output:
Array Dimension  = [5, 6, 3, 1]
Cost             = 48
Optimal Multiply = (A * (B * C))

Array Dimension  = [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Cost             = 38120
Optimal Multiply = ((((((((A * B) * C) * D) * E) * F) * G) * (H * (I * J))) * (K * L))

Array Dimension  = [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Cost             = 1773740
Optimal Multiply = (A * ((((((B * C) * D) * (((E * F) * G) * H)) * I) * J) * K))

Fortran

Translation of: Python

This is a translation of the Python iterative solution.

module optim_mod
    implicit none
contains
    subroutine optim(a)
        implicit none
        integer :: a(:), n, i, j, k
        integer, allocatable :: u(:, :)
        integer(8) :: c
        integer(8), allocatable :: v(:, :)

        n = ubound(a, 1) - 1
        allocate (u(n, n), v(n, n))
        v = huge(v)
        u(:, 1) = -1
        v(:, 1) = 0
        do j = 2, n
            do i = 1, n - j + 1
                do k = 1, j - 1
                    c = v(i, k) + v(i + k, j - k) + int(a(i), 8) * int(a(i + k), 8) * int(a(i + j), 8)
                    if (c < v(i, j)) then
                        u(i, j) = k
                        v(i, j) = c
                    end if
                end do
            end do
        end do
        write (*, "(I0,' ')", advance="no") v(1, n)
        call aux(1, n)
        print *
        deallocate (u, v)
    contains
        recursive subroutine aux(i, j)
            integer :: i, j, k

            k = u(i, j)
            if (k < 0) then
                write (*, "(I0)", advance="no") i
            else
                write (*, "('(')", advance="no")
                call aux(i, k)
                write (*, "('*')", advance="no")
                call aux(i + k, j - k)
                write (*, "(')')", advance="no")
            end if
        end subroutine
    end subroutine
end module

program matmulchain
    use optim_mod
    implicit none

    call optim([5, 6, 3, 1])
    call optim([1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2])
    call optim([1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10])
end program

Output

48 (1*(2*3))
38120 ((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))
1773740 (1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))

FreeBASIC

Translation of: VBA

This is a translation of the Python iterative solution.

Dim Shared As Integer U(), V()

Sub Aux(i As Integer, j As Integer)
    Dim As Integer k = U(i, j)
    If k < 0 Then
        Print Str(i);
    Else
        Print "(";
        Aux(i, k)
        Print "*";
        Aux(i + k, j - k)
        Print ")";
    End If
End Sub

Sub Optimize(a() As Integer)
    Dim As Integer i, j, k, c
    Dim As Integer n = Ubound(a) - 1
    Redim U(n, n), V(n, n)
    For i = 1 To n
        U(i, 1) = -1
        V(i, 1) = 0
    Next i
    
    For j = 2 To n
        For i = 1 To n - j + 1
            V(i, j) = &H7FFFFFFF
            For k = 1 To j - 1
                c = V(i, k) + V(i + k, j - k) + a(i) * a(i + k) * a(i + j)
                If c < V(i, j) Then
                    U(i, j) = k
                    V(i, j) = c
                End If
            Next k
        Next i
    Next j
    
    Print V(1, n); " ";
    Aux(1, n)
    Print
    Erase U, V
End Sub

Dim As Integer A1(1 To 4) = {5, 6, 3, 1}
Optimize(A1())
Dim As Integer A2(1 To 13) = {1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2}
Optimize(A2())
Dim As Integer A3(1 To 12) = {1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10}
Optimize(A3())

Sleep
Output:
 48(1*(2*3))
 38120((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))
 1773740(1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))

Go

The first for loop is based on the pseudo and Java code from the Wikipedia article.

package main

import "fmt"

// PrintMatrixChainOrder prints the optimal order for chain
// multiplying matrices.
// Matrix A[i] has dimensions dims[i-1]×dims[i].
func PrintMatrixChainOrder(dims []int) {
	n := len(dims) - 1
	m, s := newSquareMatrices(n)

	// m[i,j] will be minimum number of scalar multiplactions
	// needed to compute the matrix A[i]A[i+1]…A[j] = A[i…j].
	// Note, m[i,i] = zero (no cost).
	// s[i,j] will be the index of the subsequence split that
	// achieved minimal cost.
	for lenMinusOne := 1; lenMinusOne < n; lenMinusOne++ {
		for i := 0; i < n-lenMinusOne; i++ {
			j := i + lenMinusOne
			m[i][j] = -1
			for k := i; k < j; k++ {
				cost := m[i][k] + m[k+1][j] + dims[i]*dims[k+1]*dims[j+1]
				if m[i][j] < 0 || cost < m[i][j] {
					m[i][j] = cost
					s[i][j] = k
				}
			}
		}
	}

	// Format and print result.
	const MatrixNames = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
	var subprint func(int, int)
	subprint = func(i, j int) {
		if i == j {
			return
		}
		k := s[i][j]
		subprint(i, k)
		subprint(k+1, j)
		fmt.Printf("%*s -> %s × %s%*scost=%d\n",
			n, MatrixNames[i:j+1],
			MatrixNames[i:k+1],
			MatrixNames[k+1:j+1],
			n+i-j, "", m[i][j],
		)
	}
	subprint(0, n-1)
}

func newSquareMatrices(n int) (m, s [][]int) {
	// Allocates two n×n matrices as slices of slices but
	// using only one [2n][]int and one [2n²]int backing array.
	m = make([][]int, 2*n)
	m, s = m[:n:n], m[n:]
	tmp := make([]int, 2*n*n)
	for i := range m {
		m[i], tmp = tmp[:n:n], tmp[n:]
	}
	for i := range s {
		s[i], tmp = tmp[:n:n], tmp[n:]
	}
	return m, s
}

func main() {
	cases := [...][]int{
		{1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2},
		{1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10},
	}
	for _, tc := range cases {
		fmt.Println("Dimensions:", tc)
		PrintMatrixChainOrder(tc)
		fmt.Println()
	}
}
Output:
Dimensions: [1 5 25 30 100 70 2 1 100 250 1 1000 2]
          AB -> A × B           cost=125
         ABC -> AB × C          cost=875
        ABCD -> ABC × D         cost=3875
       ABCDE -> ABCD × E        cost=10875
      ABCDEF -> ABCDE × F       cost=11015
     ABCDEFG -> ABCDEF × G      cost=11017
          IJ -> I × J           cost=25000
         HIJ -> H × IJ          cost=25100
  ABCDEFGHIJ -> ABCDEFG × HIJ   cost=36118
          KL -> K × L           cost=2000
ABCDEFGHIJKL -> ABCDEFGHIJ × KL cost=38120

Dimensions: [1000 1 500 12 1 700 2500 3 2 5 14 10]
         BC -> B × C          cost=6000
        BCD -> BC × D         cost=6012
         EF -> E × F          cost=1750000
        EFG -> EF × G         cost=1757500
       EFGH -> EFG × H        cost=1757506
    BCDEFGH -> BCD × EFGH     cost=1763520
   BCDEFGHI -> BCDEFGH × I    cost=1763530
  BCDEFGHIJ -> BCDEFGHI × J   cost=1763600
 BCDEFGHIJK -> BCDEFGHIJ × K  cost=1763740
ABCDEFGHIJK -> A × BCDEFGHIJK cost=1773740

Haskell

import Data.List (elemIndex)
import Data.Char (chr, ord)
import Data.Maybe (fromJust)

mats :: [[Int]]
mats =
  [ [5, 6, 3, 1]
  , [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
  , [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
  ]

cost :: [Int] -> Int -> Int -> (Int, Int)
cost a i j
  | i < j =
    let m =
          [ fst (cost a i k) + fst (cost a (k + 1) j) +
           (a !! i) * (a !! (j + 1)) * (a !! (k + 1))
          | k <- [i .. j - 1] ]
        mm = minimum m
    in (mm, fromJust (elemIndex mm m) + i)
  | otherwise = (0, -1)

optimalOrder :: [Int] -> Int -> Int -> String
optimalOrder a i j
  | i < j =
    let c = cost a i j
    in "(" ++ optimalOrder a i (snd c) ++ optimalOrder a (snd c + 1) j ++ ")"
  | otherwise = [chr ((+ i) $ ord 'a')]

printBlock :: [Int] -> IO ()
printBlock v =
  let c = cost v 0 (length v - 2)
  in putStrLn
       ("for " ++
        show v ++
        " we have " ++
        show (fst c) ++
        " possibilities, z.B " ++ optimalOrder v 0 (length v - 2))

main :: IO ()
main = mapM_ printBlock mats
Output:
for [5,6,3,1] we have 48 possibilities, z.B (a(bc))
for [1,5,25,30,100,70,2,1,100,250,1,1000,2] we have 38120 possibilities, z.B ((((((((ab)c)d)e)f)g)(h(ij)))(kl))
for [1000,1,500,12,1,700,2500,3,2,5,14,10] we have 1773740 possibilities, z.B (a((((((bc)d)(((ef)g)h))i)j)k)

J

This is no more than a mindless transliteration of the Wikipedia Java code (for moo; for pooc, the author found Go to have the clearest expression for transliteration).

Given J's incredible strengths with arrays and matrices, the author is certain there is a much more succinct and idiomatic approach available, but hasn't spent the time understanding how the Wikipedia algorithm works, so hasn't made an attempt at a more native J solution. Others on RC are welcome and invited to do so.

moo =: verb define
  s =. m =. 0 $~ ,~ n=._1+#y
  for_lmo. 1+i.<:n do.
    for_i. i. n-lmo do.
      j =. i + lmo
      m =. _ (<i;j)} m
      for_k. i+i.j-i do.
        cost =. ((<i;k){m) + ((<(k+1);j){m) + */ y {~ i,(k+1),(j+1)
        if. cost < ((<i;j){m) do.
          m =. cost (<i;j)} m
          s =. k (<i;j)} s
        end.
      end.
    end.
  end.

  m;s
)

poco =: dyad define
  'i j' =. y
   if. i=j do.
     a. {~ 65 + i      NB. 65 = a.i.'A'
   else.
     k =. x {~ <y      NB. y = i,j
     '(' , (x poco i,k) , (x poco j ,~ 1+k) , ')'
   end.
)

optMM =: verb define
  'M S' =. moo y
  smoutput 'Cost: ' , ": x: M {~ <0;_1
  smoutput 'Order: ', S poco 0 , <:#M
)
Output:
   optMM 5 6 3 1
Cost: 48
Order: (A(BC))

   optMM 1 5 25 30 100 70 2 1 100 250 1 1000 2
Cost: 38120
Order: ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))

   optMM 1000 1 500 12 1 700 2500 3 2 5 14 10
Cost: 1773740
Order: (A((((((BC)D)(((EF)G)H))I)J)K))

Java

Thanks to the Wikipedia page for a working Java implementation.

import java.util.Arrays;

public class MatrixChainMultiplication {

    public static void main(String[] args) {
        runMatrixChainMultiplication(new int[] {5, 6, 3, 1});
        runMatrixChainMultiplication(new int[] {1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2});
        runMatrixChainMultiplication(new int[] {1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10});
    }
    
    private static void runMatrixChainMultiplication(int[] dims) {
        System.out.printf("Array Dimension  = %s%n", Arrays.toString(dims));
        System.out.printf("Cost             = %d%n", matrixChainOrder(dims));
        System.out.printf("Optimal Multiply = %s%n%n", getOptimalParenthesizations());
    }

    private static int[][]cost;
    private static int[][]order;
    
    public static int matrixChainOrder(int[] dims) {
        int n = dims.length - 1;
        cost = new int[n][n];
        order = new int[n][n];

        for (int lenMinusOne = 1 ; lenMinusOne < n ; lenMinusOne++) {
            for (int i = 0; i < n - lenMinusOne; i++) {
                int j = i + lenMinusOne;
                cost[i][j] = Integer.MAX_VALUE;
                for (int k = i; k < j; k++) {
                    int currentCost = cost[i][k] + cost[k+1][j] + dims[i]*dims[k+1]*dims[j+1];
                    if (currentCost < cost[i][j]) {
                        cost[i][j] = currentCost;
                        order[i][j] = k;
                    }
                }
            }
        }
        return cost[0][n-1];
    }

    private static String getOptimalParenthesizations() {
        return getOptimalParenthesizations(order, 0, order.length - 1);
    }
    
    private static String getOptimalParenthesizations(int[][]s, int i, int j) {
        if (i == j) {
            return String.format("%c", i+65);
        }
        else {
            StringBuilder sb = new StringBuilder();
            sb.append("(");
            sb.append(getOptimalParenthesizations(s, i, s[i][j]));
            sb.append(" * ");
            sb.append(getOptimalParenthesizations(s, s[i][j] + 1, j));
            sb.append(")");
            return sb.toString();
        }
    }

}
Output:
Array Dimension  = [5, 6, 3, 1]
Cost             = 48
Optimal Multiply = (A * (B * C))

Array Dimension  = [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Cost             = 38120
Optimal Multiply = ((((((((A * B) * C) * D) * E) * F) * G) * (H * (I * J))) * (K * L))

Array Dimension  = [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Cost             = 1773740
Optimal Multiply = (A * ((((((B * C) * D) * (((E * F) * G) * H)) * I) * J) * K))

jq

Translation of: Wren
Works with: jq

Works with gojq, the Go implementation of jq

# Input: array of dimensions
# output: {m, s}
def optimalMatrixChainOrder:
  . as $dims
  | (($dims|length) - 1) as $n
  | reduce range(1; $n) as $len ({m: [], s: []};
      reduce range(0; $n-$len) as $i (.;
        ($i + $len) as $j
        | .m[$i][$j] = infinite
        | reduce range($i; $j) as $k (.;
            ($dims[$i] * $dims [$k + 1] * $dims[$j + 1]) as $temp
            | (.m[$i][$k] + .m[$k + 1][$j] + $temp) as $cost
            | if $cost < .m[$i][$j]
              then .m[$i][$j] = $cost
              | .s[$i][$j] = $k
              else .
	      end ) )) ;

# input: {s}
def printOptimalChainOrder($i; $j):
  if $i == $j
  then [$i + 65] | implode #=> "A", "B", ...
  else "(" +
      printOptimalChainOrder($i; .s[$i][$j]) +
      printOptimalChainOrder(.s[$i][$j] + 1; $j) + ")"
  end;
 
def dimsList: [
    [5, 6, 3, 1],
    [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2],
    [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
];

dimsList[]
| "Dims  : \(.)",
  (optimalMatrixChainOrder
   | "Order : \(printOptimalChainOrder(0; .s|length - 1))",
     "Cost  : \(.m[0][.s|length - 1])\n" )
Output:
Dims  : [5,6,3,1]
Order : (AB)
Cost  : 90

Dims  : [1,5,25,30,100,70,2,1,100,250,1,1000,2]
Order : ((((((((AB)C)D)E)F)G)(H(IJ)))K)
Cost  : 37118

Dims  : [1000,1,500,12,1,700,2500,3,2,5,14,10]
Order : (A(((((BC)D)(((EF)G)H))I)J))
Cost  : 1777600

Julia

Works with: Julia version 0.6
Translation of: MATLAB

Module:

module MatrixChainMultiplications

using OffsetArrays

function optim(a)
    n = length(a) - 1
    u = fill!(OffsetArray{Int}(0:n, 0:n), 0)
    v = fill!(OffsetArray{Int}(0:n, 0:n), typemax(Int))
    u[:, 1] .= -1
    v[:, 1] .= 0
    for j in 2:n, i in 1:n-j+1, k in 1:j-1
        c = v[i, k] + v[i+k, j-k] + a[i] * a[i+k] * a[i+j]
        if c < v[i, j]
            u[i, j] = k
            v[i, j] = c
        end
    end
    return v[1, n], aux(u, 1, n)
end

function aux(u, i, j)
    k = u[i, j]
    if k < 0
        return sprint(print, i)
    else
        return sprint(print, '(', aux(u, i, k), '×', aux(u, i + k, j - k), ")")
    end
end

end  # module MatrixChainMultiplications

Main:

println(MatrixChainMultiplications.optim([1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]))
println(MatrixChainMultiplications.optim([1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]))
Output:
(38120, "((((((((1×2)×3)×4)×5)×6)×7)×(8×(9×10)))×(11×12))")
(1773740, "(1×((((((2×3)×4)×(((5×6)×7)×8))×9)×10)×11))")

Kotlin

This is based on the pseudo-code in the Wikipedia article.

// Version 1.2.31

lateinit var m: List<IntArray>
lateinit var s: List<IntArray>

fun optimalMatrixChainOrder(dims: IntArray) {
    val n = dims.size - 1
    m = List(n) { IntArray(n) }
    s = List(n) { IntArray(n) }
    for (len in 1 until n) {
        for (i in 0 until n - len) {
            val j = i + len
            m[i][j] = Int.MAX_VALUE
            for (k in i until j) {
                val temp = dims[i] * dims [k + 1] * dims[j + 1]
                val cost = m[i][k] + m[k + 1][j] + temp
                if (cost < m[i][j]) {
                    m[i][j] = cost
                    s[i][j] = k
                }
            }
        }
    }
}

fun printOptimalChainOrder(i: Int, j: Int) {
    if (i == j)
        print("${(i + 65).toChar()}")
    else {
        print("(")
        printOptimalChainOrder(i, s[i][j])
        printOptimalChainOrder(s[i][j] + 1, j)
        print(")")
    }
}

fun main(args: Array<String>) {
    val dimsList = listOf(
        intArrayOf(5, 6, 3, 1),
        intArrayOf(1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2),
        intArrayOf(1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10)
    )
    for (dims in dimsList) {
        println("Dims  : ${dims.asList()}")
        optimalMatrixChainOrder(dims)
        print("Order : ")
        printOptimalChainOrder(0, s.size - 1)
        println("\nCost  : ${m[0][s.size - 1]}\n")
    }
}
Output:
Dims  : [5, 6, 3, 1]
Order : (A(BC))
Cost  : 48

Dims  : [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order : ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost  : 38120

Dims  : [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order : (A((((((BC)D)(((EF)G)H))I)J)K))
Cost  : 1773740

Lua

-- Matrix A[i] has dimension dims[i-1] x dims[i] for i = 1..n
local function MatrixChainOrder(dims)
    local m = {}
    local s = {}
    local n = #dims - 1;
    -- m[i,j] = Minimum number of scalar multiplications (i.e., cost)
    -- needed to compute the matrix A[i]A[i+1]...A[j] = A[i..j]
    -- The cost is zero when multiplying one matrix
    for i = 1,n do
        m[i] = {}
        m[i][i] = 0
        s[i] = {}
    end

    for len = 2,n do -- Subsequence lengths
        for i = 1,(n - len + 1) do
            local j = i + len - 1
            m[i][j] = math.maxinteger
            for k = i,(j - 1) do
                local cost = m[i][k] + m[k+1][j] + dims[i]*dims[k+1]*dims[j+1];
                if (cost < m[i][j]) then
                    m[i][j] = cost;
                    s[i][j] = k; --Index of the subsequence split that achieved minimal cost
                end
            end
        end
    end
    return m,s
end

local function printOptimalChainOrder(s)
    local function find_path(start,finish)
        local chainOrder = ""
        if (start == finish) then
            chainOrder = chainOrder .."A"..start
        else
            chainOrder = chainOrder .."(" .. 
                         find_path(start,s[start][finish]) ..
                         find_path(s[start][finish]+1,finish) .. ")"
        end
        return chainOrder
    end
    print("Order : "..find_path(1,#s))
end

local dimsList = {{5, 6, 3, 1},{1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2},{1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10}}

for k,dim in ipairs(dimsList) do
    io.write("Dims  : [")
    for v=1,(#dim-1) do
        io.write(dim[v]..", ")
    end
    print(dim[#dim].."]")
    local m,s = MatrixChainOrder(dim)
    printOptimalChainOrder(s)
    print("Cost  : "..tostring(m[1][#s]).."\n")
end
Output:
Dims  : [5, 6, 3, 1]
Order : (A1(A2A3))
Cost  : 48

Dims  : [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order : ((((((((A1A2)A3)A4)A5)A6)A7)(A8(A9A10)))(A11A12))
Cost  : 38120

Dims  : [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order : (A1((((((A2A3)A4)(((A5A6)A7)A8))A9)A10)A11))
Cost  : 1773740

Mathematica / Wolfram Language

Translation of: Fortran
ClearAll[optim, aux]
optim[a_List] := Module[{u, v, n, c, r, s},
  n = Length[a] - 1;
  u = ConstantArray[0, {n, n}];
  v = ConstantArray[\[Infinity], {n, n}];
  u[[All, 1]] = -1;
  v[[All, 1]] = 0;
  Do[
   Do[
    Do[
     c = 
      v[[i, k]] + v[[i + k, j - k]] + a[[i]] a[[i + k]] a[[i + j]];
      If[c < v[[i, j]],
       u[[i, j]] = k;
       v[[i, j]] = c;
      ]
     ,
     {k, 1, j - 1}
     ]
    ,
    {i, 1, n - j + 1}
    ]
   ,
   {j, 2, n}
   ];
  r = v[[1, n]];
  s = aux[u, 1, n];
  {r, s}
  ]
aux[u_, i_, j_] := Module[{k},
  k = u[[i, j]];
  If[k < 0,
   i
   ,
   Inactive[Times][aux[u, i, k], aux[u, i + k, j - k]]
   ]
  ]
{r, s} = optim[{1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2}];
r
s
{r, s} = optim[{1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10}];
r
s
Output:
38120
(((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12)
1773740
1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11)

MATLAB

Translation of: Fortran
function [r,s] = optim(a)
    n = length(a)-1;
    u = zeros(n,n);
    v = ones(n,n)*inf;
    u(:,1) = -1;
    v(:,1) = 0;
    for j = 2:n
        for i = 1:n-j+1
            for k = 1:j-1
                c = v(i,k)+v(i+k,j-k)+a(i)*a(i+k)*a(i+j);
                if c<v(i,j)
                    u(i,j) = k;
                    v(i,j) = c;
                end
            end
        end
    end
    r = v(1,n);
    s = aux(u,1,n);
end

function s = aux(u,i,j)
    k = u(i,j);
    if k<0
        s = sprintf("%d",i);
    else
        s = sprintf("(%s*%s)",aux(u,i,k),aux(u,i+k,j-k));
    end
end
Output:
[r,s] = optim([1,5,25,30,100,70,2,1,100,250,1,1000,2])

r =

       38120


s = 

    "((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))"


[r,s] = optim([1000,1,500,12,1,700,2500,3,2,5,14,10])

r =

     1773740


s = 

    "(1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))"

Nim

Translation of: Kotlin
import sequtils

type Optimizer = object
  dims: seq[int]
  m: seq[seq[Natural]]
  s: seq[seq[Natural]]


proc initOptimizer(dims: openArray[int]): Optimizer =
  ## Create an optimizer for the given dimensions.
  Optimizer(dims: @dims)

proc findMatrixChainOrder(opt: var Optimizer) =
  ## Find the best order for matrix chain multiplication.

  let n = opt.dims.high
  opt.m = newSeqWith(n, newSeq[Natural](n))
  opt.s = newSeqWith(n, newSeq[Natural](n))

  for lg in 1..<n:
    for i in 0..<(n - lg):
      let j = i + lg
      opt.m[i][j] = Natural.high
      for k in i..<j:
        let cost = opt.m[i][k] + opt.m[k+1][j] + opt.dims[i] * opt.dims[k+1] * opt.dims[j+1]
        if cost < opt.m[i][j]:
          opt.m[i][j] = cost
          opt.s[i][j] = k


proc optimalChainOrder(opt: Optimizer; i, j: Natural): string =
  ## Return the optimal chain order as a string.
  if i == j:
    result.add chr(i + ord('A'))
  else:
    result.add '('
    result.add opt.optimalChainOrder(i, opt.s[i][j])
    result.add opt.optimalChainOrder(opt.s[i][j] + 1, j)
    result.add ')'


when isMainModule:

  const
    Dims1 = @[5, 6, 3, 1]
    Dims2 = @[1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
    Dims3 = @[1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]

  for dims in [Dims1, Dims2, Dims3]:
    var opt = initOptimizer(dims)
    opt.findMatrixChainOrder()
    echo "Dims:  ", dims
    echo "Order: ", opt.optimalChainOrder(0, dims.len - 2)
    echo "Cost:  ", opt.m[0][dims.len - 2]
    echo ""
Output:
Dims:  @[5, 6, 3, 1]
Order: (A(BC))
Cost:  48

Dims:  @[1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order: ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost:  38120

Dims:  @[1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order: (A((((((BC)D)(((EF)G)H))I)J)K))
Cost:  1773740

Perl

Translation of: Raku
use strict;
use feature 'say';

sub matrix_mult_chaining {
    my(@dimensions) = @_;
    my(@cp,@path);

    # a matrix never needs to be multiplied with itself, so it has cost 0
    $cp[$_][$_] = 0 for keys @dimensions;

    my $n = $#dimensions;
    for my $chain_length (1..$n) {
        for my $start (0 .. $n - $chain_length - 1) {
            my $end = $start + $chain_length;
            $cp[$end][$start] = 10e10;
            for my $step ($start .. $end - 1) {
                my $new_cost = $cp[$step][$start]
                             + $cp[$end][$step + 1]
                             + $dimensions[$start] * $dimensions[$step+1] * $dimensions[$end+1];
                if ($new_cost < $cp[$end][$start]) {
                    $cp[$end][$start] = $new_cost; # cost
                    $cp[$start][$end] = $step;     # path
                }
            }
       }
    }

    $cp[$n-1][0] . ' ' . find_path(0, $n-1, @cp);
}

sub find_path {
    my($start,$end,@cp) = @_;
    my $result;

    if ($start == $end) {
        $result .= 'A' . ($start + 1);
    } else {
        $result .= '(' .
                   find_path($start, $cp[$start][$end], @cp) .
                   find_path($cp[$start][$end] + 1, $end, @cp) .
                   ')';
    }
    return $result;
}

say matrix_mult_chaining(<1 5 25 30 100 70 2 1 100 250 1 1000 2>);
say matrix_mult_chaining(<1000 1 500 12 1 700 2500 3 2 5 14 10>);
Output:
38120 ((((((((A1A2)A3)A4)A5)A6)A7)(A8(A9A10)))(A11A12))
1773740 (A1((((((A2A3)A4)(((A5A6)A7)A8))A9)A10)A11))

Phix

As per the wp pseudocode

with javascript_semantics
function optimal_chain_order(int i, int j, sequence s)
    if i==j then
        return i+'A'-1
    end if
    return "("&optimal_chain_order(i,s[i,j],s)
              &optimal_chain_order(s[i,j]+1,j,s)&")"
end function
 
function optimal_matrix_chain_order(sequence dims)
    integer n = length(dims)-1
    sequence m = repeat(repeat(0,n),n),
             s = deep_copy(m)
    for len=2 to n do
        for i=1 to n-len+1 do
            integer j = i+len-1
            m[i][j] = -1
            for k=i to j-1 do
                atom cost := m[i][k] + m[k+1][j] + dims[i]*dims[k+1]*dims[j+1]
                if m[i][j]<0
                or cost<m[i][j] then
                    m[i][j] = cost;
                    s[i][j] = k;
                end if
            end for
        end for
    end for
    return {optimal_chain_order(1,n,s),m[1,n]}
end function
 
constant tests = {{5, 6, 3, 1},
                  {1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2},
                  {1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10}}
for i=1 to length(tests) do
    sequence ti = tests[i]
    printf(1,"Dims  : %s\n",{sprint(ti)})
    printf(1,"Order : %s\nCost  : %d\n",optimal_matrix_chain_order(ti))
end for
Output:
Dims  : {5,6,3,1}
Order : (A(BC))
Cost  : 48
Dims  : {1,5,25,30,100,70,2,1,100,250,1,1000,2}
Order : ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost  : 38120
Dims  : {1000,1,500,12,1,700,2500,3,2,5,14,10}
Order : (A((((((BC)D)(((EF)G)H))I)J)K))
Cost  : 1773740

Python

We will solve the task in three steps:

1) Enumerate all ways to parenthesize (using a generator to save space), and for each one compute the cost. Then simply look up the minimal cost.

2) Merge the enumeration and the cost function in a recursive cost optimizing function. The computation is roughly the same, but it's much faster as some steps are removed.

3) The recursive solution has many duplicates computations. Memoize the previous function: this yields a dynamic programming approach.

Enumeration of parenthesizations

def parens(n):
    def aux(n, k):
        if n == 1:
            yield k
        elif n == 2:
            yield [k, k + 1]
        else:
            a = []
            for i in range(1, n):
                for u in aux(i, k):
                    for v in aux(n - i, k + i):
                        yield [u, v]
    yield from aux(n, 0)

Example (in the same order as in the task description)

for u in parens(4):
    print(u)

[0, [1, [2, 3]]]
[0, [[1, 2], 3]]
[[0, 1], [2, 3]]
[[0, [1, 2]], 3]
[[[0, 1], 2], 3]

And here is the optimization step:

def optim1(a):
    def cost(k):
        if type(k) is int:
            return 0, a[k], a[k + 1]
        else:
            s1, p1, q1 = cost(k[0])
            s2, p2, q2 = cost(k[1])
            assert q1 == p2
            return s1 + s2 + p1 * q1 * q2, p1, q2
    cmin = None
    n = len(a) - 1
    for u in parens(n):
        c, p, q = cost(u)
        if cmin is None or c < cmin:
            cmin = c
            umin = u
    return cmin, umin

Recursive cost optimization

The previous function optim1 already used recursion, but only to compute the cost of a given parens configuration, whereas another function (a generator actually) provides these configurations. Here we will do both recursively in the same function, avoiding the computation of configurations altogether.

def optim2(a):
    def aux(n, k):
        if n == 1:
            p, q = a[k:k + 2]
            return 0, p, q, k
        elif n == 2:
            p, q, r = a[k:k + 3]
            return p * q * r, p, r, [k, k + 1]
        else:
            m = None
            p = a[k]
            q = a[k + n]
            for i in range(1, n):
                s1, p1, q1, u1 = aux(i, k)
                s2, p2, q2, u2 = aux(n - i, k + i)
                assert q1 == p2
                s = s1 + s2 + p1 * q1 * q2
                if m is None or s < m:
                    m = s
                    u = [u1, u2]
            return m, p, q, u
    s, p, q, u = aux(len(a) - 1, 0)
    return s, u

Memoized recursive call

The only difference between optim2 and optim3 is the @memoize decorator. Yet the algorithm is way faster with this. According to Wikipedia, the complexity falls from O(2^n) to O(n^3). This is confirmed by plotting log(time) vs log(n) for n up to 580 (this needs changing Python's recursion limit).

def memoize(f):
    h = {}
    def g(*u):
        if u in h:
            return h[u]
        else:
            r = f(*u)
            h[u] = r
            return r
    return g

def optim3(a):
    @memoize
    def aux(n, k):
        if n == 1:
            p, q = a[k:k + 2]
            return 0, p, q, k
        elif n == 2:
            p, q, r = a[k:k + 3]
            return p * q * r, p, r, [k, k + 1]
        else:
            m = None
            p = a[k]
            q = a[k + n]
            for i in range(1, n):
                s1, p1, q1, u1 = aux(i, k)
                s2, p2, q2, u2 = aux(n - i, k + i)
                assert q1 == p2
                s = s1 + s2 + p1 * q1 * q2
                if m is None or s < m:
                    m = s
                    u = [u1, u2]
            return m, p, q, u
    s, p, q, u = aux(len(a) - 1, 0)
    return s, u

Putting all together

import time

u = [[1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2],
     [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]]

for a in u:
    print(a)
    print()
    print("function     time       cost   parens  ")
    print("-" * 90)
    for f in [optim1, optim2, optim3]:
        t1 = time.clock()
        s, u = f(a)
        t2 = time.clock()
        print("%s %10.3f %10d   %s" % (f.__name__, 1000 * (t2 - t1), s, u))
    print()

Output (timings are in milliseconds)

[1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]

function     time       cost   parens
------------------------------------------------------------------------------------------
optim1    838.636      38120   [[[[[[[[0, 1], 2], 3], 4], 5], 6], [7, [8, 9]]], [10, 11]]
optim2     80.628      38120   [[[[[[[[0, 1], 2], 3], 4], 5], 6], [7, [8, 9]]], [10, 11]]
optim3      0.373      38120   [[[[[[[[0, 1], 2], 3], 4], 5], 6], [7, [8, 9]]], [10, 11]]

[1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]

function     time       cost   parens
------------------------------------------------------------------------------------------
optim1    223.186    1773740   [0, [[[[[[1, 2], 3], [[[4, 5], 6], 7]], 8], 9], 10]]
optim2     27.660    1773740   [0, [[[[[[1, 2], 3], [[[4, 5], 6], 7]], 8], 9], 10]]
optim3      0.307    1773740   [0, [[[[[[1, 2], 3], [[[4, 5], 6], 7]], 8], 9], 10]]

A mean on 1000 loops to get a better precision on the optim3, yields respectively 0.365 ms and 0.287 ms.

Iterative solution

In the previous solution, memoization is done blindly with a dictionary. However, we need to compute the optimal products for all sublists. A sublist is described by its first index and length (resp. i and j+1 in the following function), hence the set of all sublists can be described by the indices of elements in a triangular array u. We first fill the "solution" (there is no product) for sublists of length 1 (u[0]), then for each successive length we optimize using what when know about smaller sublists. Instead of keeping track of the optimal solutions, the single needed one is computed in the end.

def optim4(a):
    global u
    n = len(a) - 1
    u = [None] * n
    u[0] = [[None, 0]] * n
    for j in range(1, n):
        v = [None] * (n - j)
        for i in range(n - j):
            m = None
            for k in range(j):
                s1, c1 = u[k][i]
                s2, c2 = u[j - k - 1][i + k + 1]
                c = c1 + c2 + a[i] * a[i + k + 1] * a[i + j + 1]
                if m is None or c < m:
                    s = k
                    m = c
            v[i] = [s, m]
        u[j] = v
    def aux(i, j):
        s, c = u[j][i]
        if s is None:
            return i
        else:
            return [aux(i, s), aux(i + s + 1, j - s - 1)]
    return u[n - 1][0][1], aux(0, n - 1)


print(optim4([1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]))
print(optim4([1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]))

Output

(38120, [[[[[[[[0, 1], 2], 3], 4], 5], 6], [7, [8, 9]]], [10, 11]])
(1773740, [0, [[[[[[1, 2], 3], [[[4, 5], 6], 7]], 8], 9], 10]])

R

aux <- function(i, j, u) {
  k <- u[[i, j]]
  if (k < 0) {
    i
  } else {
    paste0("(", Recall(i, k, u), "*", Recall(i + k, j - k, u), ")")
  }
}

chain.mul <- function(a) {
  n <- length(a) - 1
  u <- matrix(0, n, n)
  v <- matrix(0, n, n)
  u[, 1] <- -1

  for (j in seq(2, n)) {
    for (i in seq(n - j + 1)) {
      v[[i, j]] <- Inf
      for (k in seq(j - 1)) {
        s <- v[[i, k]] + v[[i + k, j - k]] + a[[i]] * a[[i + k]] * a[[i + j]]
        if (s < v[[i, j]]) {
          u[[i, j]] <- k
          v[[i, j]] <- s
        }
      }
    }
  }

  list(cost = v[[1, n]], solution = aux(1, n, u))
}

chain.mul(c(5, 6, 3, 1))
# $cost
# [1] 48

# $solution
# [1] "(1*(2*3))"

chain.mul(c(1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2))
# $cost
# [1] 38120

# $solution
# [1] "((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))"

chain.mul(c(1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10))
# $cost
# [1] 1773740

# $solution
# [1] "(1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))"

Racket

Memoization

#lang racket

(define (memoize f)
  (define table (make-hash))
  (λ args (hash-ref! table args (thunk (apply f args)))))

(struct $ (cost expl))
(define @ vector-ref)

(define (+: #:combine [combine (thunk* #f)] . xs)
  ($ (apply + (map $-cost xs)) (apply combine (map $-expl xs))))

(define (min: . xs) (argmin $-cost xs))

(define (compute dims)
  (define loop
    (memoize
     (λ (left right)
       (cond
         [(= 1 (- right left)) ($ 0 left)]
         [else (for/fold ([ans ($ +inf.0 #f)]) ([mid (in-range (add1 left) right)])
                 (min: ans (+: (loop left mid) (loop mid right)
                               ($ (* (@ dims left) (@ dims mid) (@ dims right)) #f)
                               #:combine (λ (left-answer right-answer _)
                                           (list left-answer '× right-answer)))))]))))
  (loop 0 (sub1 (vector-length dims))))

Main

(define-syntax-rule (echo <x> ...)
  (begin (printf "~a: ~a\n" (~a (quote <x>) #:min-width 12) <x>) ...))

(define (solve input)
  (match-define-values ((list ($ cost explanation)) _ time _) (time-apply compute (list input)))
  (echo input time cost explanation)
  (newline))

(solve #(1 5 25 30 100 70 2 1 100 250 1 1000 2))
(solve #(1000 1 500 12 1 700 2500 3 2 5 14 10))

Output (timings are in milliseconds)

input       : #(1 5 25 30 100 70 2 1 100 250 1 1000 2)
time        : 1
cost        : 38120
explanation : ((((((((0 × 1) × 2) × 3) × 4) × 5) × 6) × (7 × (8 × 9))) × (10 × 11))

input       : #(1000 1 500 12 1 700 2500 3 2 5 14 10)
time        : 0
cost        : 1773740
explanation : (0 × ((((((1 × 2) × 3) × (((4 × 5) × 6) × 7)) × 8) × 9) × 10))

Raku

(formerly Perl 6) This example is based on Moritz Lenz's code, written for Carl Mäsak's Perl 6 Coding Contest, in 2010. Slightly simplified, it fulfills the Rosetta Code task as well.

sub matrix-mult-chaining(@dimensions) {
    my @cp;
    # @cp has a dual function: 
    # * the upper triangle of the diagonal matrix stores the cost (c) for
    #   multiplying matrices $i and $j in @cp[$j][$i], where $j > $i
    # * the lower triangle stores the path (p) that was used for the lowest cost
    #   multiplication to get from $i to $j.

    # a matrix never needs to be multiplied with itself, so it has cost 0
    @cp[$_][$_] = 0 for @dimensions.keys;
    my @path;

    my $n = @dimensions.end;
    for 1 .. $n -> $chain-length {
        for 0 .. $n - $chain-length - 1 -> $start {
            my $end = $start + $chain-length;
            @cp[$end][$start] = Inf;  # until we find a better connection
            for $start .. $end - 1 -> $step {
                my $new-cost = @cp[$step][$start]
                             + @cp[$end][$step + 1]
                             + [*] @dimensions[$start, $step+1, $end+1];
                if $new-cost < @cp[$end][$start] {
                    @cp[$end][$start] = $new-cost; # cost
                    @cp[$start][$end] = $step;     # path
                }
            }
       }
    }

    sub find-path(Int $start, Int $end) {
        if $start == $end {
            take 'A' ~ ($start + 1);
        } else {
            take '(';
            find-path($start, @cp[$start][$end]);
            find-path(@cp[$start][$end] + 1, $end);
            take ')';
        }
    }
 
   return @cp[$n-1][0], gather { find-path(0, $n - 1) }.join;
}

say matrix-mult-chaining(<1 5 25 30 100 70 2 1 100 250 1 1000 2>);
say matrix-mult-chaining(<1000 1 500 12 1 700 2500 3 2 5 14 10>);
Output:
(38120 ((((((((A1A2)A3)A4)A5)A6)A7)(A8(A9A10)))(A11A12)))
(1773740 (A1((((((A2A3)A4)(((A5A6)A7)A8))A9)A10)A11)))

Rust

use std::collections::HashMap;

fn main() {
    println!("{}\n", mcm_display(vec![5, 6, 3, 1]));
    println!(
        "{}\n",
        mcm_display(vec![1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2])
    );
    println!(
        "{}\n",
        mcm_display(vec![1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10])
    );
}

fn mcm_display(dims: Vec<i32>) -> String {
    let mut costs: HashMap<Vec<i32>, (i32, Vec<usize>)> = HashMap::new();
    let mut line = format!("Dims : {:?}\n", dims);
    let ans = mcm(dims, &mut costs);
    let mut mats = (1..=ans.1.len() + 1)
        .map(|x| x.to_string())
        .collect::<Vec<String>>();
    for i in 0..ans.1.len() {
        let mat_taken = mats[ans.1[i]].clone();
        mats.remove(ans.1[i]);
        mats[ans.1[i]] = "(".to_string() + &mat_taken + "*" + &mats[ans.1[i]] + ")";
    }
    line += &format!("Order: {}\n", mats[0]);
    line += &format!("Cost : {}", ans.0);
    line
}

fn mcm(dims: Vec<i32>, costs: &mut HashMap<Vec<i32>, (i32, Vec<usize>)>) -> (i32, Vec<usize>) {
    match costs.get(&dims) {
        Some(c) => c.clone(),
        None => {
            let ans = if dims.len() == 3 {
                (dims[0] * dims[1] * dims[2], vec![0])
            } else {
                let mut min_cost = std::i32::MAX;
                let mut min_path = Vec::new();
                for i in 1..dims.len() - 1 {
                    let taken = dims[(i - 1)..(i + 2)].to_vec();
                    let mut rest = dims[..i].to_vec();
                    rest.extend_from_slice(&dims[(i + 1)..]);
                    let a1 = mcm(taken, costs);
                    let a2 = mcm(rest, costs);
                    if a1.0 + a2.0 < min_cost {
                        min_cost = a1.0 + a2.0;
                        min_path = vec![i - 1];
                        min_path.extend_from_slice(&a2.1);
                    }
                }
                (min_cost, min_path)
            };
            costs.insert(dims, ans.clone());
            ans
        }
    }
}
Output:
Dims : [5, 6, 3, 1]
Order: (1*(2*3))
Cost : 48

Dims : [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order: ((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))
Cost : 38120

Dims : [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order: (1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))
Cost : 1773740

Stata

Recursive solution

Translation of: Python

Here is the equivalent of optim3 in Python's solution. Memoization is done with an associative array. Multiple results are returned in a structure. The same effect as optim2 can be achieved by removing the asarray machinery.

mata
struct ans {
	real scalar p,q,s
	string scalar u
}

struct ans scalar function aux(n,k) {
	external dim,opt
	struct ans scalar r,r1,r2
	real scalar s,i
	
	if (n==1) {
		r.p = dim[k]
		r.q = dim[k+1]
		r.s = 0
		r.u = strofreal(k)
		return(r)
	} else if (n==2) {
		r.p = dim[k]
		r.q = dim[k+2]
		r.s = r.p*r.q*dim[k+1]
		r.u = sprintf("(%f*%f)",k,k+1)
		return(r)
	} else if (asarray_contains(opt,(n,k))) {
		return(asarray(opt,(n,k)))
	} else {
		r.p = dim[k]
		r.q = dim[k+n]
		r.s = .
		for (i=1; i<n; i++) {
			r1 = aux(i,k)
			r2 = aux(n-i,k+i)
			s = r1.s+r2.s+r1.p*r1.q*r2.q
			if (s<r.s) {
				r.s = s
				r.u = sprintf("(%s*%s)",r1.u,r2.u)
			}
		}
		asarray(opt,(n,k),r)
		return(r)
	}
}

function optim(a) {
	external dim,opt
	struct ans scalar r
	real scalar t
	
	timer_clear()
	dim = a
	opt = asarray_create("real",2)
	timer_on(1)
	r = aux(length(a)-1,1)
	timer_off(1)
	t = timer_value(1)[1]
	printf("%10.0f %10.0f %s\n",t*1000,r.s,r.u)
}

optim((1,5,25,30,100,70,2,1,100,250,1,1000,2))
optim((1000,1,500,12,1,700,2500,3,2,5,14,10))
end

Output

         0      38120 ((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))
        16    1773740 (1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))

The timing is in milliseconds, but the time resolution is too coarse to get a usable result. A mean on 1000 loops doing the same computation yields respectively 5.772 ms and 4.430 ms for these two cases. For comparison, the computation was made on the same machine as the Python solution.

Iterative solution

Translation of: Fortran
mata
function aux(u,i,j) {
	k = u[i,j]
	if (k<0) {
		printf("%f",i)
	} else {
		printf("(")
		aux(u,i,k)
		printf("*")
		aux(u,i+k,j-k)
		printf(")")
	}
}

function optim(a) {
	n = length(a)-1
	u = J(n,n,.)
	v = J(n,n,.)
	u[.,1] = J(n,1,-1)
	v[.,1] = J(n,1,0)
	for (j=2; j<=n; j++) {
		for (i=1; i<=n-j+1; i++) {
			for (k=1; k<j; k++) {
				c = v[i,k]+v[i+k,j-k]+a[i]*a[i+k]*a[i+j]
				if (c<v[i,j]) {
					u[i,j] = k
					v[i,j] = c
				}
			}
		}
	}
	printf("%f ",v[1,n])
	aux(u,1,n)
	printf("\n")
}

optim((1,5,25,30,100,70,2,1,100,250,1,1000,2))
optim((1000,1,500,12,1,700,2500,3,2,5,14,10))
end

Output

38120 ((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))
1773740 (1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))

This solution is faster than the recursive one. The 1000 loops run now in 0.234 ms and 0.187 ms per loop on average.

VBA

Translation of: Fortran
Option Explicit
Option Base 1
Dim N As Long, U() As Long, V() As Long
Sub Optimize(A As Variant)
    Dim I As Long, J As Long, K As Long, C As Long
    N = UBound(A) - 1
    ReDim U(N, N), V(N, N)
    For I = 1 To N
        U(I, 1) = -1
        V(I, 1) = 0
    Next I
    
    For J = 2 To N
        For I = 1 To N - J + 1
            V(I, J) = &H7FFFFFFF
            For K = 1 To J - 1
                C = V(I, K) + V(I + K, J - K) + A(I) * A(I + K) * A(I + J)
                If C < V(I, J) Then
                    U(I, J) = K
                    V(I, J) = C
                End If
            Next K
        Next I
    Next J

    Debug.Print V(1, N);
    Call Aux(1, N)
    Debug.Print
    Erase U, V
End Sub
Sub Aux(I As Long, J As Long)
    Dim K As Long
    K = U(I, J)
    If K < 0 Then
        Debug.Print CStr(I);
    Else
        Debug.Print "(";
        Call Aux(I, K)
        Debug.Print "*";
        Call Aux(I + K, J - K)
        Debug.Print ")";
    End If
End Sub
Sub Test()
    Call Optimize(Array(5, 6, 3, 1))
    Call Optimize(Array(1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2))
    Call Optimize(Array(1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10))
End Sub

Output

48 (1*(2*3))
38120 ((((((((1*2)*3)*4)*5)*6)*7)*(8*(9*10)))*(11*12))
1773740 (1*((((((2*3)*4)*(((5*6)*7)*8))*9)*10)*11))

Wren

Translation of: Kotlin
var m = []
var s = []

var optimalMatrixChainOrder = Fn.new { |dims|
    var n = dims.count - 1
    m = List.filled(n, null)
    s = List.filled(n, null)
    for (i in 0...n) {
        m[i] = List.filled(n, 0)
        s[i] = List.filled(n, 0)
    }
    for (len in 1...n) {
        for (i in 0...n-len) {
            var j = i + len
            m[i][j] = 1/0
            for (k in i...j) {
                var temp = dims[i] * dims [k + 1] * dims[j + 1]
                var cost = m[i][k] + m[k + 1][j] + temp
                if (cost < m[i][j]) {
                    m[i][j] = cost
                    s[i][j] = k
                }
            }
        }
    }
}

var printOptimalChainOrder
printOptimalChainOrder = Fn.new { |i, j|
    if (i == j) {
        System.write(String.fromByte(i + 65))
    } else {
        System.write("(")
        printOptimalChainOrder.call(i, s[i][j])
        printOptimalChainOrder.call(s[i][j] + 1, j)
        System.write(")")
    }
}

var dimsList = [
    [5, 6, 3, 1],
    [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2],
    [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
]
for (dims in dimsList) {
    System.print("Dims  : %(dims)")
    optimalMatrixChainOrder.call(dims)
    System.write("Order : ")
    printOptimalChainOrder.call(0, s.count - 1)
    System.print("\nCost  : %(m[0][s.count - 1])\n")
}
Output:
Dims  : [5, 6, 3, 1]
Order : (A(BC))
Cost  : 48

Dims  : [1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2]
Order : ((((((((AB)C)D)E)F)G)(H(IJ)))(KL))
Cost  : 38120

Dims  : [1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10]
Order : (A((((((BC)D)(((EF)G)H))I)J)K))
Cost  : 1773740

zkl

Translation of: Python
fcn optim3(a){	     // list --> (int,list)
   aux:=fcn(n,k,a){  // (int,int,list) --> (int,int,int,list)
      if(n==1){
	 p,q := a[k,2];
	 return(0,p,q,k);
      }
      if(n==2){
	 p,q,r := a[k,3];
	 return(p*q*r, p, r, T(k,k+1));
      }
      m,p,q,u := Void, a[k], a[k + n], Void;
      foreach i in ([1..n-1]){
      #if 0	// 0.70 sec for both tests
	 s1,p1,q1,u1 := self.fcn(i,k,a);
	 s2,p2,q2,u2 := self.fcn(n - i, k + i, a);
      #else	// 0.33 sec for both tests
	 s1,p1,q1,u1 := memoize(self.fcn, i,k,a);
	 s2,p2,q2,u2 := memoize(self.fcn, n - i, k + i, a);
      #endif
	 _assert_(q1==p2);
	 s:=s1 + s2 + p1*q1*q2;
	 if((Void==m) or (s<m)) m,u = s,T(u1,u2);
      }
      return(m,p,q,u);
   };

   h=Dictionary();		// reset memoize
   s,_,_,u := aux(a.len() - 1, 0,a);
   return(s,u);
}

var h;		// a Dictionary, set/reset in optim3()
fcn memoize(f,n,k,a){
   key:="%d,%d".fmt(n,k);	// Lists make crappy keys
   if(r:=h.find(key)) return(r);
   r:=f(n,k,a);
   h[key]=r;
   return(r);
}
fcn pp(u){	// pretty print a list of lists
   var letters=["A".."Z"].pump(String);
   u.pump(String,
      fcn(n){ if(List.isType(n)) String("(",pp(n),")") else letters[n] })
}
fcn prnt(s,u){ "%-9,d %s\n\t-->%s\n".fmt(s,u.toString(*,*),pp(u)).println() }
s,u := optim3(T(1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2));
prnt(s,u);

s,u := optim3(T(1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10));
prnt(s,u);

optim3(T(5,6,3,1)) : prnt(_.xplode());
Output:
38,120    L(L(L(L(L(L(L(L(0,1),2),3),4),5),6),L(7,L(8,9))),L(10,11))
	-->(((((((AB)C)D)E)F)G)(H(IJ)))(KL)

1,773,740 L(0,L(L(L(L(L(L(1,2),3),L(L(L(4,5),6),7)),8),9),10))
	-->A((((((BC)D)(((EF)G)H))I)J)K)

48        L(0,L(1,2))
	-->A(BC)