converting matlab machine learning routines to php, need matrix operations

sneakyimp

I am trying to convert some matlab/octave code that implements a Support Vector Machine (SVM) into PHP. This code was distributed with homework assignments as part of Andrew Ng's 'Supervised Machine Learning: Regression and Classification' class on the Coursera website. Don't worry, I have already completed the class, and converting it to PHP is not part of the assignment.

NOTE: I've looked into FANN, but this requires the installation of another php extension. My goal here is ultimately to make lightweight, dependency-free PHP code that doesn't require exec or evocations of other languages.

For starters, I'll need to perform some matrix operations: multiply, transpose, subtract, exponentiation, etc. Here's an example matlab/octave routine. The supplied parameters x1 and x2 are row vectors:

function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp(-sum((x1 - x2) .^ 2) / (2 * sigma^2));



% =============================================================

end

I've googled 'PHP matrix operations' and the PHP libs that appear in the results are disappointing for various reasons. Is there no orthodox/canonical way of performing matrix operations in PHP?

PHP-ML is pretty impressive, and actually implements a Support Vector Classifier. However, installing it brings in 150 PHP files, and the SVM is implemented by calling exec on an executable. Installing gives you linux, windows, and macos executables.

MarkBaker's PHPMatrix kinda works but its constructor seems to convert row vectors to column vectors and vice versa. This:

$matrix = new Matrix\Matrix([1,2,3,4]);
print_r($matrix);

appears to construct a column vector:

Matrix\Matrix Object
(
    [rows:protected] => 4
    [columns:protected] => 1
    [grid:protected] => Array
        (
            [0] => Array
                (
                    [0] => 1
                )

            [1] => Array
                (
                    [0] => 2
                )

            [2] => Array
                (
                    [0] => 3
                )

            [3] => Array
                (
                    [0] => 4
                )

        )

)

This code creates the exact same column vector, with 4 rows and one column:

$matrix = new Matrix\Matrix([[1],[2],[3],[4]]);

It gets even more confusing because matrix multiplication is not a symmetric operation. Consider this matlab code:

m1 = [1,1,1,1]
m2 = [1,3,5,7]

m1' * m2
ans =
   1   3   5   7
   1   3   5   7
   1   3   5   7
   1   3   5   7

m2 * m1'
ans =  16

You'd think this code using the MarkBaker lib would yield that first 4x4 matrix, but it does not:

$m1 = new Matrix\Matrix([1,1,1,1]);
$m2 = new Matrix\Matrix([1,3,5,7]);
$v = $m1->transpose()->multiply($m2);
print_r($v->toArray());

The result:

Array
(
    [0] => Array
        (
            [0] => 16
        )

)

While one might rationalize the order of operations to say that $m1 is being applied as a multiplication operation to $m2, it's just confusing.

And then there's NumPHP. The documentation for extremely sparse and not especially helpful. No details are given about transposition, for instance. I also can't seem to find any way to perform matrix multiplication -- the documentation mentions only the dot product. Unless I'm mistaken, matrix multiplication and dot products are not the same thing.

Any suggestions for a canonical/orthodox/intuitive approach to matrix operations would be much appreciated.

Weedpacket

sneakyimp Unless I'm mistaken, matrix multiplication and dot products are not the same thing.

I had a look at the code and dot does a matrix multiplication.

Basically, ordinary vector dot product is a special case of matrix multiplication. (Snipping a bunch of tensor algebra that boils down to:...) treating the first vector a 1×n matrix (a "covector") and the second as n×1 (a "vector"), their product is a 1×1 scalar result.

It's an area where there is a lot of old and confusing notation that is still hanging around and still being taught.

sneakyimp

Weedpacket I had a look at the code and dot does a matrix multiplication.

Could you clarify 'does a matrix multiplication?'

In octave/matlab, this won't work, and elicits a complaint:

x1 = [1,1,1,1]
x2 = [1,3,5,7]
x1 * x2
error: operator *: nonconformant arguments (op1 is 1x4, op2 is 1x4)

Whereas numphp has no problem doing this:

$m1 = new NumArray(
        [1,1,1,1]
);
echo $m1 . "\n";
$m2 = new NumArray(
        [1,3,5,7]
);
echo $m2 . "\n";
$v = $m1->dot($m2);
echo "$v\n";

The output:

NumArray([1, 1, 1, 1])
NumArray([1, 3, 5, 7])
NumArray(16)

It seems like dot($somevar) will transpose $somevar automatically when doing the multiplication, regardless of whether that is what is intended or not. This might be handy in the majority of cases (I would not know) but it seems problematic if you don't want to do a dot product and are actually trying to do a strict matrix multiplication.

I apologize if I seem obtuse here, but I am truly confused. I would again point out that matrix multiplication is not a symmetric operation:

x1 = [1,1,1,1];
x2 = [1,3,5,7];

x1 * x2'
ans =  16

x2' * x1
ans =

   1   1   1   1
   3   3   3   3
   5   5   5   5
   7   7   7   7

I don't see how you can get those different results using NumPHP:

$m1 = new NumArray([1,1,1,1]);
$m2 = new NumArray([1,3,5,7]);
$v = $m1->dot($m2);
echo "$v\n";
// NumArray(16)


$m1 = new NumArray([1,1,1,1]);
$m2 = new NumArray([1,3,5,7]);
$v = $m2->dot($m1);
echo "$v\n";
// NumArray(16)

$m1 = new NumArray([1,1,1,1]);
$m2 = new NumArray([1,3,5,7]);
$v = $m1->dot($m2->getTranspose());
echo "$v\n";
// NumArray(16)

$m1 = new NumArray([1,1,1,1]);
$m2 = new NumArray([1,3,5,7]);
$v = $m2->dot($m1->getTranspose());
echo "$v\n";
// NumArray(16)

$m1 = new NumArray([1,1,1,1]);
echo "$m1\n";
// NumArray([1, 1, 1, 1])
$m2 = $m1->getTranspose();
echo "$m1\n";
// NumArray([1, 1, 1, 1])

Given how confusing matrix operations can be, it seems helpful to have your code complain when you try to multiply mismatched matrices. It's also helpful to know if your matrix is a column vector or a row rector. It seems very fishy that calling getTranspose() returns something identical to the matrix you called it on.

EDIT: still more confusing:

$m1 = new NumArray([1,1,1,1]);
echo "$m1\n";
$m2 =new NumArray([[1],[3],[5],[7]]);
echo "$m2\n";


$v = $m1->dot($m2);
echo "$v\n";

output:

NumArray([1, 1, 1, 1])

NumArray([
  [1],
  [3],
  [5],
  [7]
])

NumArray([16])

Weedpacket

Well... maybe I shouldn't have snipped the stuff I snipped.

Thing is, every tensor has a rank, scalars are rank 0, vectors have rank 1, matrices have rank 2, and so on; basically it's a matter of how deeply the array brackets are nested. Dot product/matrix multiplication mushes the last dimension of its first argument with the first dimension of its second argument. [1,1,1,1] isn't a 1×n matrix because it only has rank 1; a 1×n matrix has rank two (with dimensions 1 and n) and would be [[1,1,1,1]]. Nor is it a (rank-6) 1×1×1×1×1×n tensor. Or 1×1×1×n×1×1 for that matter.

Mushed together: elements are multiplied pairwise and the resulting products summed. For a rank-n tensor T and a rank-m tensor U their product is:

T                 · U                 = ∑ T              U
 t1,t2,t3,...,tn     u1,u2,u3,...,um    k  t1,t2,t3,...k  k,u2,u3,...,um

Where there's one of those equations for each possible assignment of array indices to t1, t2, …, um, and k ranges over array indices as well. Which means the product of a rank n tensor and a rank m tensor is a tensor with rank n+m-2 (indexed by t1,t2,...t(n-1),u2,u3...um; in a sense one rank goes away when you stop double-counting the dimension along which the two tensors were aligned and elements paired off, and the other gets lost in the summation).

Since two vectors are both rank 1, their product has rank 1+1-2=0 — it's just a number. Since two matrices are both rank 2, their product has rank 2+2-2=2. And the product of a vector and a matrix has rank 1+2-2=1; i.e. a vector.

And, of course, implicit in the definition of tensor product is the fact that tn has to range over the same set of array indices as u1, so that k is well-defined.

Transpose? It (typically) swaps the last two dimensions of the tensor (though any permutation counts as a "transpose"). Which means for it to have any effect the thing being transposed has to have rank at least 2. "Typically" because it's usually only seen applied to matrices which only have two dimensions to swap.

So ... try multiplying [[1,1,1,1]] by [[1],[3],[5],[7]] and vice versa and compare the results with the matlab ones (does it use different operators for matrix multiplication and vector dot product? What does it do with tensors of higher rank?). Note that, as the bracket nesting shows, you're multiplying matrices, not vectors. (Incidentally, in your tests just a little further up, you echoed $m1 instead of $m2 in the last test. I would be interested to know what PHPNum reckons the transpose of a vector looks like.)

sneakyimp

I appreciate you taking the time for the detailed explanation. It has been a very long time since I took a class that bothered to discuss tensors or matrix operations in any detail, but I get the impression that this mathematical concept involves a level of desirable abstraction that goes a step beyond my prior exposure to matrix operations. What I need at the moment is for my PHP code to replicate the behavior I see in my octave/matlab code. It seems to me the distinction between vectors and matrices doesn't really exist in octave/matlab. I.e., octave treats [1,1,1,1] as a 1x4 matrix:

x1 = [1,1,1,1]
size(x1)
ans =

   1   4

I have always had great difficulty performing these operations in my head and juggling transpositions and arithmetic operations and their resulting dimensions is nearly impossible for me -- I wonder if perhaps it might reveal some undiagnosed dyslexia on my part.

Weedpacket So ... try multiplying [[1,1,1,1]] by [[1],[3],[5],[7]] and vice versa and compare the results with the matlab ones

This PHP code:

$m1 = new NumArray([[1,1,1,1]]);
echo "$m1\n";
$m2 =new NumArray([[1],[3],[5],[7]]);
echo "$m2\n";
$v = $m1->dot($m2);
echo "$v\n";

Yields this output:

NumArray([
  [1, 1, 1, 1]
])

NumArray([
  [1],
  [3],
  [5],
  [7]
])

NumArray([
  [16]
])

Octave barfs on the results, as it should, because these matrices are both 1x4:

x1 = [[1,1,1,1]]
x2 = [[1],[3],[5],[7]]
x1 * x2
error: operator *: nonconformant arguments (op1 is 1x4, op2 is 1x4)

This barfing is a helpful guide rail for me, and helps me to understand that I need a transposition, or have my dimensions reversed. On the other hand, when such an operation succeeds, this tends to suggest that I am multiplying the right numbers together.

Octave/matlab does in fact use the asterisk for matrix multiplication and has a dot function to perform a dot product. It is noteworthy that the documentation mentions vectors here:

dot (x, y, dim)
Compute the dot product of two vectors.
If x and y are matrices, calculate the dot products along the first non-singleton dimension.
If the optional argument dim is given, calculate the dot products along this dimension.
Implementation Note: This is equivalent to sum (conj (X) .* Y, dim), but avoids forming a temporary array and is faster. When X and Y are column vectors, the result is equivalent to X' * Y. Although, dot is defined for integer arrays, the output may differ from the expected result due to the limited range of integer objects.

Weedpacket Incidentally, in your tests just a little further up, you echoed $m1 instead of $m2 in the last test. I would be interested to know what PHPNum reckons the transpose of a vector looks like.

Sorry for the mistake. I can confirm this weirdness. This PHP code:

$m1 = new NumArray([1,3,5,7]);
echo "$m1\n";
$m2 = $m1->getTranspose();
echo "$m2\n";

Yields this result:

NumArray([1, 3, 5, 7])

NumArray([1, 3, 5, 7])

It appears that the two matrices are internally identical. This PHP code:

$m1 = new NumArray([1,3,5,7]);
print_r($m1);
$m2 = $m1->getTranspose();
print_r($m2);

Shows two identical objects:

NumPHP\Core\NumArray Object
(
    [shape:protected] => Array
        (
            [0] => 4
        )

    [data:protected] => Array
        (
            [0] => 1
            [1] => 3
            [2] => 5
            [3] => 7
        )

    [cache:protected] => Array
        (
        )

)
NumPHP\Core\NumArray Object
(
    [shape:protected] => Array
        (
            [0] => 4
        )

    [data:protected] => Array
        (
            [0] => 1
            [1] => 3
            [2] => 5
            [3] => 7
        )

    [cache:protected] => Array
        (
        )

)

Weedpacket

sneakyimp Octave barfs on the results, as it should, because these matrices are both 1x4:
x1 = [[1,1,1,1]]
x2 = [[1],[3],[5],[7]]
x1 * x2
error: operator *: nonconformant arguments (op1 is 1x4, op2 is 1x4)

Interesting: because it looks like x1 is 1x4 (it has one element that contains four elements) but x2 is 4x1 (it has four elements that contain one element each). So octave's just ignoring that? Ah:

If x and y are matrices, calculate the dot products along the first non-singleton dimension.

Yes, yes it is. "A two dimensional matrix where one dimension has only one index is treated as though it were a vector." Makes me wonder how it distinguishes row and column vectors from each other). But it does mean that this distinction between dot product versus matrix multiplication suddenly has to be a thing. How do transposed vectors get represented as an octave literal? It apparently happens because x2' works where x2 doesn't.

sneakyimp Sorry for the mistake. I can confirm this weirdness. This PHP code:

Makes sense (or, at least, is consistent). It's a vector; it only has one dimension, while the transpose needs at least two dimensions to work with (so at least a matrix), seeing as it's expected to swap them.

Incidentally, I looked at the NumPY code that NumPHP tries to mimic; apparently there the tensor product contracts the last dimension of the left argument with the second-to-last dimension of the right argument (instead of the first). Probably for performance's sake. For matrices it wouldn't make a difference (since the second-to-last dimension out of two is the first), but for anything of higher rank you'd be doing a bunch of before-and-after transposing to get the conventional behaviour (so much for performance). I don't know if NumPHP does that as well, but I guess it only matters if you're working with tensors of rank higher than 2.

sneakyimp

It occurred to me that perhaps some of my confusion here may be related to notation -- we use square brackets in PHP to denote arrays but in octave we use them to denote matrices or vectors. Code that does one thing in PHP might do something entirely different in octave. For example, this results in some 1x3 matrices being concatenated into a 1x9 matrix:

x1 = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
x1 =

   1   2   3   4   5   6   7   8   9

in octave you delimit rows with semicolons:

x1 = [[1, 2, 3]; [4, 5, 6]; [7, 8, 9]]
x1 =

   1   2   3
   4   5   6
   7   8   9

And the treatment of brackets, commas, and semicolons suggests to me that there's no concept of any matrix/tensor/object with a rank greater than 2:

x1 = [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]]
x1 =

   1
   2
   3
   4
   5
   6
   7
   8
   9

Weedpacket Makes sense (or, at least, is consistent). It's a vector; it only has one dimension, while the transpose needs at least two dimensions to work with (so at least a matrix), seeing as it's expected to swap them.

I have a dim inkling that this is perhaps too clever by half to say that there's no such thing as a matrix with a singleton dimension. Or perhaps there is a distinction between a vector with 4 elements and a 1x4 matrix or 4x1 matrix. Octave certainly seem to make the distinction, and complains when dimensions don't match.

Regarding performance, I have some very disappointing news for these PHP libs. I have a data object distilled from a corpus of ham & spam messages. The training set, Xtrain, has 1375 elements, which of which represents one email message in my corpus. Each element is an array with 1966 elements, each of which is a zero or one indicating whether a particular word (from a vocabulary of 1966 words) is present in that particular email message. I munged my data corpus with some other scripts and exported this 2-D array into a json file, which can be loaded thusly:

$data_file = '../../../machine-learning/training_data_sets.json';
$data = json_decode(file_get_contents($data_file), TRUE);
$keys = array_keys($data); // Xtrain, ytrain, Xval, yval, Xtest, ytest
echo "Xtrain has " . sizeof($data['Xtrain']) . " elements\n"; // 1375 elements
echo "Xtrain[0] has " . sizeof($data['Xtrain'][0]) . " elements\n"; // 1966 elements

One of the most time-consuming calculations in the svm training script is to take this Xtrain corpus and multiply it by its transpose. In octave, this happens almost instantly with an older, slightly different corpus:

size(X)
ans =

   1382   1996

% this returns insantly
K = X * X';

I've tried this using all 3 of the matrix libs I've stumbled across, and they run MUCH slower

NumPHP:

$X = new NumArray($data['Xtrain']);
$start = microtime(TRUE);
$XT = $X->getTranspose();
echo "getTranspose in " . (microtime(TRUE) - $start) . " seconds\n"; // 5.9611051082611 seconds
$start = microtime(TRUE);
$K = $X->dot($XT);
echo "dot product in " . (microtime(TRUE) - $start) . " seconds\n"; // 466.75154590607 seconds

MarkBaker:

$X = new Matrix\Matrix($data['Xtrain']);
$start = microtime(TRUE);
$XT = $X->transpose();
echo "transpose in " . (microtime(TRUE) - $start) . " seconds\n"; // 0.13952708244324 seconds
$start = microtime(TRUE);
$K = $X->multiply($XT);
echo "multiply in " . (microtime(TRUE) - $start) . " seconds\n"; // 1821.0363881588 seconds

PHP-ML:

$X = new Matrix($data['Xtrain']);
$start = microtime(TRUE);
$XT = $X->transpose();
echo "transpose in " . (microtime(TRUE) - $start) . " seconds\n"; // 0.14473080635071 seconds
$start = microtime(TRUE);
$K = $X->multiply($XT);
echo "multiply in " . (microtime(TRUE) - $start) . " seconds\n"; // 292.54629206657 seconds

I've yet to investigate whether these short PHP scripts are calculating the correct multiplication, but the performance seems really slow, especially given the speed at which octave seems to perform the calculation. Notational confusion aside, I am now questioning whether these libs are suitable at all for this application.

Weedpacket

sneakyimp I have a dim inkling that this is perhaps too clever by half to say that there's no such thing as a matrix with a singleton dimension. Or perhaps there is a distinction between a vector with 4 elements and a 1x4 matrix or 4x1 matrix. Octave certainly seem to make the distinction, and complains when dimensions don't match.

Well, it's doing it by using different punctuation to indicate different levels; none of your octave samples to now mentioned that. So [1,3,5,7], [[1,3,5,7]] and [[1],[3],[5],[7]] are all the same thing to octave and its matrix multiplication operation treats them all as 1×4 matrices; since they're all the same object as var as octave is concerned, dot is treating them all the same.

So dot([1,1,1,1],[1,3,5,7]) == [1,1,1,1]*[1;3;5;7], right? Will that latter give you a number or a 1x1 matrix? Or, for that matter, a length-1 vector? How would you even distinguish them? Either way, I expect the numerical value would be the same. What do dot([1;1;1;1],[1;3;5;7]), dot([1,1,1,1],[1;3;5;7]), and dot([1;1;1;1],[1,3,5,7]) produce?

(Incidentally, personally I've been using a different notation again; I've been translating it into PHP and what I thought was octave before learning about ; for the sake of not adding another notation into the mix.)

Basically, since octave doesn't do tensors and its 1-row matrices look exactly like vectors (and its 1-column matrices look exactly like column vectors), it's ended up with two distinct functions to do the same job (one treats [1,3,5,7] as a vector and one as a matrix); to get the same results, one of them requires you to distinguish "row" vectors and "column" vectors by using , vs. ;, and the other doesn't bother to distinguish "row" vectors and "column" vectors even if you do.

The NumPHP and NumPY libraries don't do this confounding-vectors-with-matrices thing; they retain rank information. As intended, their dot product function works on tensors with arbitrary rank — including vectors and matrices. The dot product of two vectors is a scalar (rank 1+1-2) and the dot product of a 1×n matrix with an n×1 matrix is a 1x1 matrix (rank 2+2-2) with its sole element having the same numerical value as the scalar. The same function also of course handles multiplying length-n vectors by n×m matrices and n×m by m, with the result in both cases being a vector.

sneakyimp

Weedpacket Makes sense (or, at least, is consistent). It's a vector; it only has one dimension, while the transpose needs at least two dimensions to work with (so at least a matrix), seeing as it's expected to swap them.

I have a dim inkling that this is perhaps too clever by half to say that there's no such thing as a matrix with a singleton dimension.

Yup. Octave's dot function chooses to pretend they don't exist; its matrix-multiplication operator doesn't. (The other problem is that by ignoring ranks higher than 2, it's stuck with having to come up with some other mechanism to distinguish "row vectors" and "column vectors".)

sneakyimp Octave certainly seem to make the distinction

Only if you use , and ; as appropriate to distinguish between "row" and "column" vectors (or 1×n and n×1 matrices, which to octave are the same thing). Otherwise it doesn't make the distinction which is why it needs two different functions to do the same job.

sneakyimp Regarding performance, I have some very disappointing news for these PHP libs.... Octave (a FOSS variant of matlab) handles matrix operations effortlessly.

Well, octave was built with maths processing in mind, and its matrix handling is probably heavily optimised. PHP ... isn't. (Actually, that splitting of dot product and matrix multiplication, and limiting itself to rank-2, may have been deliberate performance decisions.)

Even NumPY resorts to compiled machine code when it can (and the GPU when it can find one). If you were to go the compiled route you'd probably want to make (or find) a library that you can use PHP's FFI extension with.

NogDog

sneakyimp but the performance seems really slow

Time to switch to Rust? https://youtu.be/oY0XwMOSzq4
😉

I haven't dived into it yet, but a guy on our team is doing amazing things with it to manage the processing of huge JSON files of health insurance data.

sneakyimp

NogDog
Ugh. I confess that rust sounds interesting, and I quite enjoy hearing from you, NogDog, but that narrator claims to have written backends in python, ruby, and node, then digresses into a rant about the linux kernel adoption being proof of rust's awesomeness. Then digresses further to discuss Jack@twitter, Discord, and the perils of garbage collection and high-level languages. His accent, which exudes a kindly-but-pedantic hauteur, is infuriating. He mentions haskell, julia, and scala, but not once are matrix operations mentioned.

Octave (a FOSS variant of matlab) handles matrix operations effortlessly. I might offload the matrix operations to octave via exec, but I might as well rely on the PHP-ML lib if I must offload my calculations to some other programming language. My stated goal is to make a dependency-free PHP implementation for use in plugins for a variety of popular frameworks/CMSes. I might further refine that goal to avoid use of exec-type functions.

sneakyimp

Weedpacket So [1,3,5,7], [[1,3,5,7]] and [[1],[3],[5],[7]] are all the same thing to octave and its matrix multiplication operation treats them all as 1×4 matrices

Yes.

Weedpacket So dot([1,1,1,1],[1,3,5,7]) == [1,1,1,1]*[1;3;5;7], right?

Yes. Both of those appear to return a scalar value. HOWEVER, if you reverse the order of operands, they return different results:

octave:82> dot([1,1,1,1],[1,3,5,7])
ans = 16
octave:83> [1,1,1,1]*[1;3;5;7]
ans = 16
octave:84> dot([1,3,5,7], [1,1,1,1])
ans = 16
octave:85> dot([1;3;5;7], [1,1,1,1])
ans = 16
octave:86> [1;3;5;7]*[1,1,1,1]
ans =

   1   1   1   1
   3   3   3   3
   5   5   5   5
   7   7   7   7

Weedpacket What do dot([1;1;1;1],[1;3;5;7]), dot([1,1,1,1],[1;3;5;7]), and dot([1;1;1;1],[1,3,5,7]) produce?

octave:86> dot([1;1;1;1],[1;3;5;7])
ans = 16
octave:87> dot([1,1,1,1],[1;3;5;7])
ans = 16
octave:88> dot([1;1;1;1],[1,3,5,7])
ans = 16

Weedpacket (Incidentally, personally I've been using a different notation again; I've been translating it into PHP and what I thought was octave before learning about ; for the sake of not adding another notation into the mix.)

I appreciate your efforts, and have only slowly come to realize where my confusion lies. You are enormously helpful, and hardly to blame for Octave's quirks or dissimilarities to PHP--or my confused & anemic descriptions.

Weedpacket ended up with two distinct functions to do the same job

The dot function seems like a convenience function more than anything. I haven't seen it used in the code I've been looking at. Like many PHP functions, it surely has its own quirks to which one must become accustomed.

Weedpacket The NumPHP and NumPY libraries don't do this confounding-vectors-with-matrices thing; they retain rank information

This may indeed be the desired behavior from a sort of abstract mathematical perspective, but I would point out that matrix multiplication is not symmetric -- see those first dot-versus-multiplication results at the top of this post. I apologize that I'm unable to foresee the broader mathematical/code significance of this. My comprehension of matrix operations and their relation to the underlying algegbra is dim, at best. I remain confused about what it means that dot always returns a scalar (as NumPY presumably does) when multiplying these vectors/singleton matrices, whereas the octave/matlab multiplication operation returns a scalar in one case and a 4x4 matrix when the operands are reversed.

Weedpacket Well, octave was built with maths processing in mind, and its matrix handling is probably heavily optimised. PHP ... isn't. (Actually, that splitting of dot product and matrix multiplication, and limiting itself to rank-2, may have been deliberate performance decisions.)

Quite frankly, I am far more interested in performance for the task at hand which has two main parts:

convert this svm code from octave/matlab PHP such that it doesn't require installation of additional PHP or exec functions
insure that performance is reasonable for common use cases

I'm far less interested in any mathematical orthodoxy as far as distinctions between a vector and a 1xn or nx1 matrix. In fact, you might say I'm sort of accustomed to the quirky octave behaviors and rather hooked on the stubborn refusal of the octave multiplier operation to proceed if the dimensions don't match. When the code makes assumptions about what I'm trying to multiply, this, to me, introduces ambiguity. I like that [1,1,1,1] * [1;3;5;7] and [1;3;5;7]*[1,1,1,1] produce different results. When they don't, I feel like things are getting unpredictable.

There can be no question that octave is highly optimized for multiplying a matrix times its transpose -- at least compared to these PHP libs I've identified. Octave is three or four orders of magnitude faster on a 1375x1966 matrix. I wrote my own multiply function, and it takes over ten minutes on this laptop to multiply this matrix times its transpose. The transposition happens in about 0.14 seconds. It's the multiplication of 1375x1966 matrix by 1966x1375 matrix that takes so long. My function:

        public function multiply($matrix) {
                $iteration_count = $this->colCount();
                if ($iteration_count !== $matrix->rowCount()) {
                        throw new Exception('nonconformant arguments (op1 is ' . $this->sizeAsString() . ', op2 is ' . $matrix->sizeAsString() . ')');
                }

                $result_rows = $this->rowCount();
                $result_cols = $matrix->colCount();

                $ret_array = [];
                for($i=0; $i<$result_rows; $i++) {
                        $ret_array[$i] = array_fill(0, $result_cols, 0);
                        for($j=0; $j<$result_cols; $j++) {
                                for($k=0; $k<$iteration_count; $k++) {
                                        $ret_array[$i][$j] += ($this->data[$i][$k] * $matrix->data[$k][$j]);
                                }
                        }
                }

                return new static($ret_array);
        }

multiply in 734.63768601418 seconds

I looked into the PHP-ML code, which had the fastest multiplication operation in my earlier testing, and appropriated its multiplication algorithm to produce this function:

        public function multiply2($matrix) {
                // this fn largely borrowed from php-ml, which said:
                /*
                - To speed-up multiplication, we need to avoid use of array index operator [ ] as much as possible( See #255 for details)
                - A combination of "foreach" and "array_column" works much faster then accessing the array via index operator
                */

                $iteration_count = $this->colCount();
                if ($iteration_count !== $matrix->rowCount()) {
                        throw new Exception('nonconformant arguments (op1 is ' . $this->sizeAsString() . ', op2 is ' . $matrix->sizeAsString() . ')');
                }

                $result_cols = $matrix->colCount();
                $product = [];
                foreach ($this->data as $row => $rowData) {
                        for ($col = 0; $col < $result_cols; ++$col) {
                                $columnData = array_column($matrix->data, $col);
                                $sum = 0;
                                foreach ($rowData as $key => $valueData) {
                                        $sum += $valueData * $columnData[$key];
                                }
                                $product[$row][$col] = $sum;
                        }
                }

                return new self($product);
        }

It's nearly 3 times faster:

multiply2 in 288.33355784416 seconds

However, as I mentioned previously, octave performs a calculation like this instantly. If I had to guess, I'd say octave is probably capitalizing on some kind of contiguous memory structure where there's no need to traverse/search an array to locate a particular index. I expect there might also be some shrewd mathematical optimizations -- there's probably a symmetry that results when you multiply a matrix by its transpose, and perhaps they're using some kind of middle-out optimization or something. I seriously doubt it'll be possible to massage these PHP multiplication functions to boost performance by 3 orders of magnitude, but I'm open to suggestions. This is a drag.

Weedpacket

sneakyimp but I would point out that matrix multiplication is not symmetric -- see those first dot-versus-multiplication results at the top of this post.

This isn't really an issue if you keep in mind that vectors are not matrices; the same formula applies for vectors and matrices.
For vectors it reduces to

T    · U   = ∑ T  U
 t1     u1   k  k  k

And U·T will produce the same result

For matrix×matrix it's

 T      · U      = ∑ T     U
  t1,t2    u1,u2   k  t1,k  k,u2

which isn't the same thing as

 U      · T      = ∑ U     T
  u1,u2    t1,t2   k  u1,k  k,t2

because you're pairing different elements together (in the first, the elements in one row of T with the elements in a column of U, and in the second the elements in a row of U with the elements in a column of T). That's why matrix multiplication isn't symmetric, and also why it doesn't matter for vector dot product; vectors don't have that many dimensions so there is less freedom for how their elements can be paired up.

sneakyimp I like that [1,1,1,1] * [1;3;5;7] and [1;3;5;7]*[1,1,1,1] produce different results. When they don't, I feel like things are getting unpredictable.

Which is what you'd expect from multiplying a 1×4 matrix and a 4×1 matrix vs. a 4×1 matrix and a 1×4 matrix. Not having that ,/; distinction I'd write them as

{{1,1,1,1}}.{{1},{3},{5},{7}} = {{16}}

and

{{1},{3},{5},{7}}.{{1,1,1,1}} = {{1, 1, 1, 1}, {3, 3, 3, 3}, {5, 5, 5, 5}, {7, 7, 7, 7}}

And the brackets are significant because without them they're not matrices but vectors, which means

{1, 1, 1, 1} . {1, 3, 5, 7}={1, 3, 5, 7} . {1, 1, 1, 1}=16

sneakyimp . I remain confused about what it means that dot always returns a scalar

That I think comes back to dot deliberately ignoring singleton dimensions; it throws away the distinction between [1,1,1,1] and [1;1;1;1] which means dot([1;1;1;1].[1,1,1,1]) and dot([1,1,1,1].[1;1;1;1]) are indistinguishable; so no wonder they produce the same answer. The matrix multiplication operator doesn't do this.
Octave's notation makes it impossible to distinguish between a length-n vector and a 1×n matrix, so it has to have different functions where one treats it as a vector and the other treats it as a matrix. The one that treats it as a vector has to also treat n×1 matrices as vectors as well.

sneakyimp

Weedpacket (in the first, the elements in one row of T with the elements in a column of U, and in the second the elements in a row of U with the elements in a column of T). That's why matrix multiplication isn't symmetric, and also why it doesn't matter for vector dot product; vectors don't have that many dimensions so there is less freedom for how their elements can be paired up.

Thanks once again, Weedpacket, for patiently explaining these things. I do have a much better grip, I think, on vector/matrix operations. I feel stronger in my conceptual understanding of what is transpiring in these operations. My only lingering question is perhaps should we always assume 1xn or nx1 matrices are vectors, or are there circumstances where we'd want to treat them as actual matrices? All the code and explanations I've seen here suggest that 1xn/nx1 matrices are always intended as vectors.

And, lastly, I think it seems clear that these three PHP libs we've examined are sorely inadequate for a pure PHP implementation (no FANN, no exec of other executable or script) of my Support Vector Machine. My spam filter example requires one to multiple a 1375x1966 matrix by its inverse 1966x1375 matrix and it is a hundred times too slow. Unless someone can suggest a way to construct a data structure in pure PHP that can efficiently access matrix indices, PHP is going to be doing some searching/traversal of these large arrays. I suspect (but am not sure) that the PHP slowness here is due to array search/traversal, whereas in C or Octave (or perhaps Rust), one might be able to construct an efficient, contiguous memory segment where accessing $matrix[$i][$j] is simply a matter of calculating a memory offset from some starting memory location and then it's simple RAM access rather than array search/traversal.

Weedpacket

sneakyimp is perhaps should we always assume 1xn or nx1 matrices are vectors, or are there circumstances where we'd want to treat them as actual matrices? All the code and explanations I've seen here suggest that 1xn/nx1 matrices are always intended as vectors.

Er ... if you mean my code and explanations then that's the opposite of what I was saying.

sneakyimp My spam filter example requires one to multiple a 1375x1966 matrix by its inverse 1966x1375 matrix

You mean "transpose", of course, because multiplying a matrix by its inverse is really easy if you assume it has one 🙂

I do have a couple of suggestions; one is to use an SPLFixedArray (especially since the array's sizes are all known).
The other, if AA' is the main slow part, is to use an algorithm optimised for that case. Mainly, since the result is symmetric, you only need to compute half of it; the other half would be a mirror image. You could either calculate the half and then copy everything into the other half (if you're using SPLFixedArrays this won't make much difference memory-wise, but the time spent on it may save more time later on), or just leave the array half-filled and redirect lookups as appropriate (so, every time you want to look up $symm[$j][$i] you'd check to see if $i>$j and if it is you'd look up $symm[$i][$j] instead.)

sneakyimp

Weedpacket You mean "transpose", of course, because multiplying a matrix by its inverse is really easy if you assume it has one 🙂

OOF yes I meant transpose. Inverse is something else entirely. Sorry for the mistake.

I will of course attempt a multiply function using SPLFixedArray.

sneakyimp

Not having much luck wtih SplFixedArray. Here's a script:

class JMatrix {
	public $data;

	public function __construct($arr) {
		if (is_a($arr, 'SplFixedArray')) {
			// shortcut...provide SplFixedArray param, bypass checks and just assign
			$this->data = $arr;
		} elseif (!is_array($arr)) {
			// if not SplFixedArray, must be standard array
			throw new Exception('Constructor param must be an array');
		} else {
			// if it's an array, check if 2D array (matrix) or just 1D array vector, which we treat as matrix
			if (!is_array($arr[0])) {
				// the supplied parameter is a simple array, so we'll construct a row vector
				// from its elements
				$row_count = 1;
				$col_count = sizeof($arr);
				$this->data = new SplFixedArray(1);
				$this->data->offsetSet(0, new SplFixedArray($col_count));
				foreach($arr as $j => $v) {
					$this->data->offsetGet(0)->offsetSet($j, $v);
				}

			} else {
				// an array of arrays...sanity check on element counts
				$row_count = sizeof($arr);
				$col_count = sizeof($arr[0]);
				foreach($arr as $v) {
					if (sizeof($v) !== $col_count) {
						throw new Exception('Invalid matrix array. Inconsistent column count in constructor array.');
					}
				}
				$this->data = new SplFixedArray($row_count);
				foreach($arr as $i => $row) {
					$this->data->offsetSet($i, new SplFixedArray($col_count));
					foreach($row as $j => $v) {
						$this->data->offsetGet($i)->offsetSet($j, $v);
					}
				}
				
			}
		}
	}

	public function rowCount() {
		return $this->data->getSize();
	}
	public function colCount() {
		return $this->data->offsetGet(0)->getSize();
	}

	public function size() {
		return [$this->rowCount(), $this->colCount()];
	}
	public function sizeAsString() {
		return $this->rowCount() . 'x' . $this->colCount();
	}

	public function __toString() {
		$retval = "[\n";
		foreach($this->data as $i => $row) {
			$retval .= "\t[" . implode(",", (array)$row) . "]\n";
		}
		$retval .= "]\n";
		return $retval;
	}

	public function multiply($matrix) {
		$iteration_count = $this->colCount();
		if ($iteration_count !== $matrix->rowCount()) {
			throw new Exception('nonconformant arguments (op1 is ' . $this->sizeAsString() . ', op2 is ' . $matrix->sizeAsString() . ')');
		}

		$result_rows = $this->rowCount();
		$result_cols = $matrix->colCount();

		// accumulate results in this array
		$ret_array = new SplFixedArray($result_rows);

		// temp vars to avoid object property lookup
		$td = $this->data;
		$md = $matrix->data;

		for($i=0; $i<$result_rows; $i++) {
			// current row from this matrix
			$drow = $td->offsetGet($i);

			// accumulate output matrix row here
			$retrow = new SplFixedArray($result_cols);

			// loop thru $matrix columns
			for($j=0; $j<$result_cols; $j++) {

				// accumulate product sum here
				$sum = 0;

				// iteration count is number of columns in this matrix == number of rows in $matrix
				for($k=0; $k<$iteration_count; $k++) {
					$sum += ($drow->offsetGet($k) * $matrix->data->offsetGet($k)->offsetGet($j));
				}

				// all products summed, store it in our output matrix row
				$retrow->offsetSet($j, $sum);
			}

			// put the output matrix row into our output matrix
			$ret_array->offsetSet($i, $retrow);
		}
		return new static($ret_array);
	}

	public function multiply2($matrix) {
		$iteration_count = $this->colCount();
		if ($iteration_count !== $matrix->rowCount()) {
			throw new Exception('nonconformant arguments (op1 is ' . $this->sizeAsString() . ', op2 is ' . $matrix->sizeAsString() . ')');
		}

		$result_rows = $this->rowCount();
		$result_cols = $matrix->colCount();

		// accumulate results in this array
		$ret_array = new SplFixedArray($result_rows);

		// temp vars to avoid object property lookup
		$td = $this->data;
		$md = $matrix->data;


		$td->rewind();
		$i = 0;
		while ($td->valid()) {
			// current row from this matrix
			$drow = $td->current();

			// accumulate output matrix row here
			$retrow = new SplFixedArray($result_cols);

			// loop thru $matrix columns
			for($j=0; $j<$result_cols; $j++) {

				// accumulate product sum here
				$sum = 0;

				// iteration count is number of columns in this matrix == number of rows in $matrix
				$md->rewind();
				$k = 0;
				while($md->valid()) {
					$mrow = $md->current();
					$sum += ($drow->offsetGet($k) * $mrow->offsetGet($j));

					$md->next();
					$k++;
				}

				// all products summed, store it in our output matrix row
				$retrow->offsetSet($j, $sum);
			}

			// put the output matrix row into our output matrix
			$ret_array->offsetSet($i, $retrow);

			$td->next();
			$i++;
		}
		return new static($ret_array);
	}

	public function transpose() {
		$arr = [];
		foreach($this->data as $row) {
			$arr[] = (array)$row;
		}
		$arr = array_map(null, ...$arr);
		return new static($arr);
	}
}


echo "SANITY CHECK\n";
$a = new JMatrix([
	[1, 4, -2],
	[3, 5, -6]
]);
echo $a;
echo "size: " .  $a->sizeAsString() . "\n";

$b = new JMatrix([
	[5, 2, 8, -1],
	[3, 6, 4, 5],
	[-2, 9, 7, -3]
]);
echo $b;
echo "size: " .  $b->sizeAsString() . "\n";

$v = $a->multiply2($b);
echo $v . "\n";
echo "size: " .  $v->sizeAsString() . "\n";

$vt = $v->transpose();
echo "$vt\n";
echo "size: " .  $vt->sizeAsString() . "\n";


$data_file = './training_data_sets.json';
$data = json_decode(file_get_contents($data_file), TRUE);

echo 'Xtrain rows: ' . sizeof($data['Xtrain']) . "\n";
echo 'Xtrain cols: ' . sizeof($data['Xtrain'][0]) . "\n";

$X = new JMatrix($data['Xtrain']);

$start = microtime(TRUE);
$XT = $X->transpose();
echo "transpose in " . (microtime(TRUE) - $start) . " seconds\n";

$start = microtime(TRUE);
$K = $X->multiply2($XT);
echo "multiply2 in " . (microtime(TRUE) - $start) . " seconds\n";

multiply, my first attempt, runs in 973 seconds:

Xtrain rows: 1375
Xtrain cols: 1966
transpose in 0.79622912406921 seconds
multiply in 973.16395616531 seconds

multiply2, where i attempt to avoid indexed lookups, actually takes longer -- about 20 minutes 😢.

Xtrain rows: 1375
Xtrain cols: 1966
transpose in 0.79264497756958 seconds
multiply2 in 1199.6731460094 seconds

Weedpacket

Oh, well. Worth a try. (Incidentally, the instanceof operator could be used instead of the is_a function. Since SplFixedArray implements Countable you can use count() instead of calling ->getSize() manually; implementing IteratorAggregate means that the ->valid() ->rewind(), ->current(), and ->next() methods will be called by foreach as and when needed; and implementing the ArrayAccess interface means the ability to use array index notation as implemented by ->offsetGet() and ->offsetSet(), e.g. $this->data->offsetGet(0)->offsetSet($j, $v); becomes $this->data[0][$j] = $v;, And I'm guessing you're using PHP 7: SplFixedArray dropped the rewind() method in v8 when it was switched from being an Iterator to an IteratorAggregate. Woo. Big sidebar. But I did get some improvement by making those changes. Possibly because of the dedicated opcodes for those things.)

But I think PHP will (even with JIT enabled) lag behind domain-specific languages that will have had dedicated native code for those tasks for which the language is intended ("a matrix multiplication? I have machine code for that compiled from C++ right here"); PHP's execution would still be passing through the various reference-counted zval structures and hashtables and whatnot before getting right to the gist of the task. (See, e.g. Internal value representation in PHP 7 - Part 2.)

sneakyimp

Weedpacket Oh, well. Worth a try. (Incidentally, the instanceof operator could be used instead of the is_a function...I did get some improvement by making those changes. Possibly because of the dedicated opcodes for those things.)

Thanks for these helpful clarifications. It had occurred to me that the transpose operation is quite fast, and part of the slowness seems to be due to the array indexing required. Note how the fastest PHP implementation, PHP-ML, makes use of the array_column function and explicitly tries to avoid array indexing as much as possible. It therefore occurred to me that I might quickly transpose the matrix supplied to my multiply function and benefit from slightly more efficient array operations. However, I don't think these little changes will bring us anywhere near the supremely fast performance I see in Octave. I exported my 1375x1966 array to an octave-friendly text format and loaded it up on octave and multipled it with this little octave/matlab script:

Xtrain = dlmread("Xtrain.txt");
sz = size(Xtrain)
t0 = clock();
K = Xtrain * Xtrain';
elapsed_time = etime(clock(), t0);
fprintf("transpose multiplication completed in %d seconds\n", elapsed_time);

Octave completes the matrix multiplication in a stunning 67ms:

sz =

   1375   1966

transpose multiplication completed in 0.0668411 seconds

This is 4376.7 times faster than the fastest PHP code I've found, which is accomplished by the matrix classes in the PHP-ML library.

NOTE that the PHP-ML library apparently punts to some kind of compiled executable to perform its SVM training operations. It constructs a command and invokes the external executable with exec.

I consulted with a friend who's a talented python jockey who offered some interesting initial observations on matrix multiplication algorithms:

Strassen's was the first notable speedup, ₁₉₇₀ — O(n ^ 2.8) iirc — and people have knocked the exponent down further and further in the ensuing decades. In fact just within the last two months I read about this recent improvement: as per https://en.wikipedia.org/wiki/Matrix_multiplication#Computational_complexity,

As of December 2020, the best matrix multiplication algorithm is by Josh Alman and
Virginia Vassilevska Williams and has complexity O(n ^ 2.3728596).

That statement just quoted has a footnote to this paper: https://arxiv.org/abs/2010.05846, "A Refined Laser Method and Faster Matrix Multiplication".
It's a tiny speedup over the previous best-known algo: best was formerly

O(n^{2.37287 + ε}) for all ε;

the paper improves that bound to

O(n^{2.37286 + ε}) for all ε;

My friend rightly observes:

It would not be easy to whip up working code from the paper, written for actual and aspiring specialists.

He pointed out an interesting exercise in matrix operation optimization and noted that numpy is very fast. Like the paper, these articles get into a lot of the nitty gritty (i.e., the weeds) of bottlenecks and various optimizations.

I guess Octave and Matlab use BLAS (Basic Linear Algebra Subprograms), originally a Fortran library now ported to C or something, under the hood.

My kind friend imported my 1375x1966 matrix for some python profiling. Some interesting results on his 2015 Macbook Pro:

Numpy takes 2.06s to matmult Xtrain by its transpose

He took the trouble to write a "pure python" function which he ran in a "notebook" :

PS I'm trying to time the naive algorithm in pure Python, operating either on the same numpy "ndarray"s, or on native Python data structures (list of lists of numbers)
I've tried both, a few times; I give up waiting, I haven't seen either complete, it takes minutes (at least!) and I get bored.

Safe conclusion: Performance of pure Python naive algorithm in a notebook is worse than PHP. Way worse: in a practical sense it doesn't even terminate.
Performance at command line: (minutes later, still waiting… I'll check back in an hour, if it terminates )

Maybe it's faster at the command line than in a notebook? I bet not enough to matter, but let's try.
With native Python datatypes: it takes at least tens of seconds... a minute… 5+ minutes…
Just checked back: it took 20m 29.15s
(So, yes: faster than when running in a notebook, but still unusable.)

He went on to make some minor improvements, but the 'pure python' approach, which doesn't make use of NumPY, is clearly slow like the pure PHP approach. Oddly enough, I find this something of a relief. I'd hate to think python is inherently better than my dear old PHP.

The lesson here is that it's clearly advantageous--necessary!--to delegate matrix operations to something outside PHP (or python) which takes advantage of decades of algorithmic refinement by the world's finest nerds. I truly doubt that your typical linux server will have octave installed, but it's easy enough with apt if you have su access via CLI. I think your typical server is much more likely to have python installed. Or perhaps I could distribute my own executable to handle matrix operations? Any suggestions about farming out these matrix operations from PHP would be much appreciated. Ideally we'd find some tool for matrix operations that will already exist in your typical server so we can just exec it from PHP with reasonable hope of success.

Weedpacket

The last time I had to do anything numerically-matricy (it was actually a dynamic programming problem) that was just way too slow and/or large to run in PHP, I whipped up a program in C to do it and compiled (TinyCC was good enough) it as an executable that I communicated with via proc_open. That was years and major versions ago. If PHP had had FFI back then I would have considered compiling it as a library and calling the API directly (warning: requires third-party library, which I suspect is the same one NumPY ends up using) instead of that mucking around with stdio. (The PHP manual has a note that accessing FFI structures is twice as slow as accessing PHP ones, so "it makes no sense to use the FFI extension for speed". But accessing the structure isn't the bottleneck here.)

sneakyimp

Thanks, Weedpacket, for the detail. Your post is over my head, and I'll have to do some remedial study on the particulars of compiling/linking/assembling to understand the dynamics here. Offhand, it seems tricky to find some speedy lib external to PHP which can perform fast matrix operations. I expect there are two paths:
1) try using exec to run some CLI commands in the hope that python or octave or some other package has been installed on the machine. Might also have to download a python lib or two like NumPy.
2) try invoking some executable file that I cook up myself via exec or FFI or some other process-to-process method. This seems quite involved, and it sounds like a real chore to distribute executables that'll work on linux, MacOS, windows, or whatever.

Obviously, my goal of a 'pure PHP' library is unattainable. I still hope that I might work up a solution that doesn't require su access to install (e.g., to configure PHP with FFI module or FANN, or to install octave or python).