Matrix by vector multiplication

[Michael Bell]

Simple...?

Matrix by vector multiplication using MAT:

Say, for example, M() is a matrix of 3 rows by 3 columns, V() is a column vector of 3 rows...

DIM M(1 TO 3,1 TO 3) AS EXT, V(1 TO 3) AS EXT

If, after populating these arrays with data, I wish to replace the data in V() with the result of multiplying M() by V(), as in...

MAT V() = M() * V()

I get an incorrect result, whereas...

MAT V2() = M() * V()

provides the correct result, but I then have to copy V2() into V().

Providing a second dimension in V(), as in V(1 to 3,1 to 1), doesn't help.

What am I missing here? I'm perfectly prepared for someone to tell me that I'm being very stupid!

Michael Bell (using PB/CC5)

[Arie Verheul]

Matrix multiplication

Matrix multiplication is very different from ordinary multiplication.
It is not so easy to explain in a few lines, but maybe you should see a text book on this subject.

Later: there is something else with PB matrices :

All operations involving two or more arrays require that they be of exactly the same size and type, without exception. Failure to adhere causes undefined results.

Arie Verheul

[Aldo Vitagliano]

Let the 3 component of the vector V be X, Y, Z, and the elements of the matrix M be M11, M12, M13, M21, M22, M23, M31, M32, M33 (1st index is the row).

Then the three component of the producy vector are:
M11*X + M12*Y + M13*Z
M21*X + M22*Y + M23*Z
M31*X + M32*Y + M33*Z

Most likely, PB computes one component at a time, storing it in the array declared to contain the result. If this is the same as the starting V array, when the second component is calculated, the wrong value of X is used, and so on ...

[Michael Bell]

Matrix by vector multiplication

Thanks Aldo and Arie.

I don't have a problem with matrix multiplication as such, and in fact find it very useful to couch various statistical and mathematical models in compact matrix form. My problem is, I think, that I am not properly understanding where PB is putting the intermediate values in the calculation before returning the result. For the statement MAT V()=M()*V(), I was expecting the result vector to be formed in some temporary space before being copied over the original vector, whereas - I'm sure you are right here Aldo - it looks like the intermediate quantities are copied into the original vector in the midst of the overall operation (although, intuitively, this should at least lead to the first vector element being correct, whereas even this is wrong...).

Hence, MAT V()=M()*V() yields an incorrect result, whereas MAT V2()=M()*V() performs as expected.

The issue seems not to be one of matrix size. I think that HELP is somewhat misleading here - what matters is that the matrices should be conformable rather than identical in size, hence it is entirely admissible to multiply a 3x3 matrix with a column vector of 3 values. Just to check this out, I tried padding out the column vector into a 3x3 matrix with zeroes on the last two columns. I get exactly the same results as before - the same incorrect values if the results matrix variable also appears on the right-hand side of the calculation, but the correct values if I use a new matrix to hold the results.

So... this leads me to treat MAT statements very carefully, and I may be inclined to use my own (slower) functions rather than MAT statements when I am uncertain about what PB is doing behind the scenes. Should the hazards of using the same array on both sides of the calculation in MAT statements be made clearer in HELP? Or else is this something that should be remedied...? Using, for example, the R statistical programming language, a statement such as V<-M%*%V, where M is a matrix and V is a compatible vector (<- does the assignment and %*% does the matrix multiplication) works just fine.

[Chris Holbrook]

Should the hazards of using the same array on both sides of the calculation in MAT statements be made clearer in HELP? Or else is this something that should be remedied...?

My guess is that Powerbasic Inc Support could answer those questions, but you will have to ask them direct.

[Michael Bell]

Thanks Chris. I have contacted PowerBASIC Support. The answer is that a third array is always needed for matrix multiplication using MAT statements. I have suggested that this be made clear in HELP, since at present the only relevant statement is "it is permissible to specify one array for multiple MAT parameter", which sort of implies that it would be OK to use the same array on both sides of an expression. One would naturally assume that coding x=y*x would be fine for scalar arithmetic, and it's not going too far to assume the same for matrix arithmetic, as is true for many high-level statistical and mathematical packages. Updating a vector by repeated multiplication by a matrix is a common operation in some types of numeric application, so if a statement is available to perform such operations it is important to be clear about how it operates.

[Chris Holbrook]

The old Help system is getting a bit of a hammering at the moment, isn't it?

If I understand you correctly then you need to use a temporary matrix variable when doing this multiplication - I wonder if you could declare it as a macrotemp variable and do the multiplication with a macro? That way it would feel more like a native function than if you code it out every time.

[Aldo Vitagliano]

I have an additional observation about Matrix products.
If speed is an issue, then it looks better to directly write the multiplication code in the program than to use the bult-in MAT routine.
Here is a sample code comparing the two different approaches:

[FONT="Courier New"]
DEFDBL A-G, M, O-Z
DEFLNG H-L, N
FUNCTION PBMAIN()
   REGISTER J AS LONG, K AS LONG
   N=1000
   DIM M(N,N), Q(N,N), Q1(N,N), Q2(N,N)
   RANDOMIZE TIMER
   '... Matrix fill-up
   FOR I=1 TO N: FOR J=1 TO N: M(I,J)=RND: Q(I,J)=I+J: NEXT: NEXT
   T0=TIMER
   '... Multiply M()*Q()
   FOR I=1 TO N
      FOR J=1 TO N
        SUM=0
        FOR K=1 TO N: SUM=SUM+M(I,K)*Q(K,J): NEXT
        Q1(I,J)=SUM
      NEXT J
   NEXT I
   T1=TIMER-T0: T0=TIMER
   '... Built-in multiplication
   MAT Q2()=M()*Q()
   T2=TIMER-T0
   PRINT "Coded product    Built-in product"
   PRINT USING$ ("####.#s            ", T1,T2)
   PRINT
   '... Correctness test
   OK&=-1
   FOR I=1 TO N: FOR J=1 TO N: IF Q2(I,J)<>Q1(I,J) THEN OK&=0
   NEXT: NEXT
   IF OK& THEN PRINT "OK!" ELSE PRINT "Comparison failed!"
   WAITKEY$
END FUNCTION

[/FONT]

[Michael Bell]

Chris
Thanks, good suggestion, and it will force me to learn to use macros.

Aldo
Very interesting - directly written code seems to be more than twice as fast as the MAT statement. Copying the results matrix back over the original makes negligible difference to the overall times, whether a MAT statement or directly written element-by-element copying is used.

[Robert DeBolt]

Matrix Inverse calculation

I ran Aldo's program on my 366mhz intel thinkpad and got the following timings:

'Matrix Fill-up time =  .691000000000582
'
'Coded product    Built-IN product
' 127.3s             286.6s
'
' 286.6 / 127.3 = 2.25
'
'OK!

'Matrix Fill-up time =  .721000000000386
'
'Coded product    Built-IN product
' 139.9s             295.8s
'
' 295.8 / 139.9 = 2.11
'
'OK!

That's Impressive!
So I thought about benchmarking the Matrix Inverse function, which is probably the most calculation intensive function, and got the following result:

'
'Coded inverse    Built-IN inverse
' 360.4s            1038.7s
'
' 1038.7 / 360.45 = 2.88
'
'OK!
'

Now that is impressive!!
I translated the FreeBasic function "gausjord" from the "fbmath" library.
Google on "fbmath" and download from the http://sourceforge.net/projects/fbmath/ site.

The source codes and exec are in the attached zip file.