Pascal - Matrix Chain Multiplication

[Muttley]

Hello, i have some problems implementing the following pseudo-code. Te problem is that in Pascal i can return user defined data structures (in this case: matrix), and i also tried to return it as a parameter with "var" but it gave a stack overflow error. Anyone can help me?

pseudo (the well known matrix chain multiplication algorithm):

MCM (X,s,i,j)
  if j>i then m1=MCM(X,s,i,s[i,j])
                 m2=MCM(X,s,s[i,j]+1,j)
                 return MatMult(m1,m2)
          else return X(i)

X - sequence of matrixes
s - split table
MatMult - multiplies 2 matrixes (i wrote it so that it stores the result in a third one like MatMult(m1,m2,m3)

my pascal implementation:

procedure MCM(x1:seqmat; s1:matrix; i1,j1:integer; var res:matrix);
  var a,b,c:matrix;
  begin
    if j1>i1 then begin
	               MCM(x1,s1,i1,s1.matr[i1,j1],a);
		       MCM(x1,s1,s1.matr[i1,j1+1],j1,b);
		       MatMult(a,b,res);
                       end
		else res:=x1[i1];
  end;

[Duoas]

OK, since no one else is looking at this one... (I had to go learn to do the Matrix Chain Multiplication)

I think you are missing something here.
The MCM algorithm isn't actually supposed to multiply the matrices together.
It is just to tell you the best order in which to multiply the matrices.

Once you know the order, then you can do the multiplications.

The other trick is a sneaky one: that sl matrix you've got there is a dirty trick to replace a binary tree. Except it is missing something: links to the next node. Meaning you need another matrix to track that information.

To review the algorithm, then, using our matrix/binary-tree substitutes, it should work something like the following example illustrates.

Given: A1 x A2 x A3 x A4
where

• A1 is 2 rows by 3 columns
• A2 is 5 rows by 10 columns
• A3 is 10 rows by 7 columns
• A4 is 50 rows by 5 columns

Question: What is the best way to parenthesize this equation? (Or, which two matrices do I multiply together first, then which one do I multiply with that, etc)

(Sorry, but I don't know how to do subscripts with DaniWeb, so pretend that the 1,2,3,4 are subscripts to A.)

Since we've got four matrices to order, we need two 4 by 4 matrices: one to represent the best cost of multiplying some set, and another to indicate where to go next. Let's build that now:

least costs                   next node
   1    2    3    4              1    2    3    4
1 ____ ____ ____ ____         1 ____ ____ ____ ____
2 ____ ____ ____ ____         2 ____ ____ ____ ____
3 ____ ____ ____ ____         3 ____ ____ ____ ____
4 ____ ____ ____ ____         4 ____ ____ ____ ____

Now, it doesn't actually matter which way you order this, but I will say that the column represents the starting subscript and the row represents the ending subscript. We can immediately cross the diagonal and the top half of each matrix out --we won't need those positions. Also, the leaf nodes don't have next nodes, so we can ignore those too:

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 ____                        2
3 ____ ____                   3 ____
4 ____ ____ ____              4 ____ ____

Alright, let's start figuring out best order.

First we'll take them by pairs:
A1 x A2 costs 2 x 3 x 10 = 60 multiplications. Let's note that in our table. column 1 to row 2:

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 _[U]60[/U]_                        2
3 ____ ____                   3 ____
4 ____ ____ ____              4 ____ ____

A2 x A3 costs 5 x 10 x 7 = 450 multiplications. Let's add that: column 2 to row 3:

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 _[U]60[/U]_                        2
3 ____ _[U]450[/U]                   3 ____
4 ____ ____ ____              4 ____ ____

And finally, A3 x A4 costs 10 x 7 x 5 = 450. Column 3 to row 4:

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 _[U]60[/U]_                        2
3 ____ _[U]450[/U]                   3 ____
4 ____ ____ _[U]450[/U]              4 ____ ____

Now it is time to fill-in the costs for larger subsequences. Before we took them by pairs. Now we'll take them by triplets.
A1 x A2 x A3 costs either (A1 x A2) x A3 or A1 x (A2 x A3). Looking in our table we can fill in some of our information, and get: (col1,row2) x A3 or A1 x (col2,row3), which becomes either (60) x 7 = 420 or 3 x (450) = 1350. Choosing the cheapest of these means picking the first.

So (A1 x A2) x A3 costs 420 multiplications. Let's add that to our tables, and make a note of which matrix subscript is outside the parentheses:

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 _[U]60[/U]_                        2
3 _[U]420[/U] _[U]450[/U]                   3 __[U]3[/U]_
4 ____ ____ _[U]450[/U]              4 ____ ____

Doing the same for A2 x A3 x A4, we find that (A2 x A3) x A4 = 2250 and A2 x (A3 x A4) = 4500. This time the cheapest is with A4 on the outside:

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 _[U]60[/U]_                        2
3 _[U]420[/U] _[U]450[/U]                   3 __[U]3[/U]_
4 ____ [U]2250[/U] _[U]450[/U]              4 ____ __[U]4[/U]_

Finally, we can calculate the total cheapest cost for the whole thing (taking it four at a time):
(A1 x A2 x A3) x A4 = (420) x 5 = 2100
A1 x (A2 x A3 x A4) = 3 x (2250) = 6750
The clear winner is the first.

least costs                   next node
   1    2    3    4              1    2    3    4
1                             1
2 _[U]60[/U]_                        2
3 _[U]420[/U] _[U]450[/U]                   3 __[U]3[/U]_
4 [U]2100[/U] [U]2250[/U] _[U]450[/U]              4 __[U]4[/U]_ __[U]4[/U]_

These two things are the return from the MCM algorithm:

The cheapest cost is 2100 multiplications needed to multiply the three matrices, and the order, outside-in is: A4, A3, (A1 A2)

Now, hopefully that helps with MCM. Now to give quick consideration to doing it recursively. Because in the recursive method, we start at the root of the tree and work our way out to the leaves, no? Whereas above we went from leaves to root.

Well, this is the power of recursion. Starting at the root, we want to know which has the better cost: A1 x (A2 x A3 x A4) or (A1 x A2 x A3) x A4. Since one of the two multiplicands is compound, we need to recurse, twice. Once to find the MCM for (A2 x A3 x A4) and once for (A1 x A2 x A3), which we can then compare and choose against.

We stop recursing when both terms of the equation represent single matrices. So (A1 x A2) is simply (rows of A1) x (cols of A1) x (cols of A2), with no previous. We stick that result in our charts and return it without recursing.

The power is that, for each subset (or node in our tree), we only need to calculate the total cost once, since we can refer to the tree to determine if we have already found a value. So, in the above example, we only needed to calculate the total cost of A2 x A3 once, instead of twice.

Whew. Hope this helps.

[Muttley]

Ok i get the point, but there are two possibilities:
1. i didn't understood the whole matrix chain multiplication
2. or u didn't understood my question

my ideea :

We have a sequence of matrices and we need to multiply them in an optimla way.
The first thing to do is to store the dimension of the matrix in an array (p) like that:
A(i) is p(i-1) x p(i)
ex: A(1) 3x5 A(2) 5x8 A(3) 8x6 A(4) 6x2
p(0)=3 p(1)=5 p(2)=8 p(3)=6 p(4)=2 (length(p)=number of matrices+1)

After that we construct two matrices m1 and s1 according to the sequence of dimensions,
where m1 is the minimum cost table, and s1 is the split table (so i think this gives the order)

procedure MCO (p1:sequence; var s1,m1:matrix);
  var l,r,k,i,j:byte; q:integer;
  begin
    m1.rows:=n;
    m1.columns:=n;
    s1.rows:=n;
    s1.columns:=n;
    for i:=1 to n do
	m1.matr[i,i]:=0;
    for l:=2 to n do
      for i:=1 to n-l+1 do
	  begin
          j:=i+l-1;
          m1.matr[i,j]:=MAXINT;
          for k:=i to j-1 do
		begin
              q:=m1.matr[i,k]+m1.matr[k+1,j]+p1[i-1]*p1[k]*p1[j];
              if (q<m1.matr[i,j]) then begin
				   		     m1.matr[i,j]:=q;
						     s1.matr[i,j]:=k;
                                       end;
            end;
        end;
  end;

explication: matrix - user defined data type containing the matrix itself and its dimensions

according to the split table (s1) now we have the order of the multiplication.
printig the order (optimal parens) on the screen:

procedure POP(s1:matrix; i1,j1:byte);
  begin
    if i1=j1 then write('A',i1)
             else begin
                    write('(');
                    POP(s1,i1,s1.matr[i1,j1]);
                    POP(s1,s1.matr[i1,j1]+1,j1);
                    write(')');
			end;
  end;

we just have to implement the subalgorithm that using the split table will multiply the matrices, which in pseudo code looks like this:

Matrix-Chain-Multiply(A,s1,i,j)
if j > i
then
x = Matrix-Chain-Multiply(A,s1,i,s1[i,j]);
y = Matrix-Chain-Multiply(A,s1,s1[i,j]+1,j);
return MatMult(x,y);
else
return A(i)

explication: A - the sequence of matrices
MatMult - simply multiplies two matrices and returns the result

the question is how to implement in Pascal this pseudo code, because the one that i give above gives stack overflow error.

[Duoas]

2. It is entirely possible that I have misunderstood something about your question. It isn't actually very clear. I know you are frustrated by a stack overflow with your posted algorithm. I have largely ignored that part because:

1. You have definitely misunderstood the MCM.

MCM does not give you the quotient matrix. It only computes the optimal way to multiply the matrices together.

So your MCO procedure corresponds to the actual Matrix Chain Multiplication algorithm. Except that it doesn't properly work.

Firstly, p is incorrect. You cannot ignore the number of rows for each matrix but the first. Also, your loops look a little wonky to me. I haven't traced through them, but there are a couple of things that make them very suspect to me. I suggest you spend some time reconsidering how they work.

Secondly, your last two procs assume something about the data in the split table -- it acts properly on only one of the two possibilities.

The reason I went through all that above is because I want you to forget the pseudo-code stuff someone else gave you and to think about the way the algorithm works on your own. I want you to do this because it is what you need to do to really fix the underlying problems you are having with the code.

As to the stack overflow, if your matrices are of any considerable dimensions then each one chews up a very large amount of space on the stack. You can fix it by making a matrix a specific object (which gets allocated on the heap, not the stack) or by using an out parameter (a var parameter used for output).

Hope this helps.

[Muttley]

The MCM must multiply and give back the result of multiplying the sequence of matrices ... it's a recursive subalgorithm that calls itselg until the data that it returns will be a matrix of the sequence ... and we also call the MatMult algorithm that multiplies two matrices ... otherwise it doesnt makes sence ....

Another thing p i correct, we can ignore the rows of the matrices because two matrices can be multiplied if and only if the number of columns of the first matrix equals the number of rows of the second matrix

[Duoas]

It makes perfect sense. You just need to think about it a little more. I don't intend to post back again just to argue what MCM is supposed to do.

As for p, I understand how matrix multiplication works. Don't be fooled by the nonsense numbers I used above in the example. I suppose you can do it that way, but I think you are making more work for yourself than you have to...

Hope this helps.