Added convolve operation to matrix and vector

Author: andrewdalpino

CHANGELOG.md

@@ -1,6 +1,7 @@
 - Unreleased
     - Added ref using row elimination method
     - Added universal comparison methods to tensor API
+    - Added convolve operation to Vector and Matrix
 
 - 1.0.1
     - Added Column Vector

src/Matrix.php

@@ -771,18 +771,22 @@ public function matmul(Matrix $b) : self
 
         $bT = $b->transpose();
 
-        $c = [[]];
+        $c = [];
+
+        foreach ($this->a as $row) {
+            $temp = [];
 
-        foreach ($this->a as $i => $row) {
             foreach ($bT as $column) {
                 $sigma = 0;
 
-                foreach ($row as $k => $value) {
-                    $sigma += $value * $column[$k];
+                foreach ($row as $i => $value) {
+                    $sigma += $value * $column[$i];
                 }
 
-                $c[$i][] = $sigma;
+                $temp[] = $sigma;
             }
+
+            $c[] = $temp;
         }
 
         return self::quick($c);
@@ -798,13 +802,62 @@ public function matmul(Matrix $b) : self
     public function dot(Vector $b) : ColumnVector
     {
         if ($this->n !== $b->size()) {
-            throw new DimensionalityMismatchException("Matrix A requires"
+            throw new DimensionalityMismatchException('Matrix A requires'
                 . " $this->n elements but Vector B has {$b->size()}.");
         }
 
         return $this->matmul($b->asColumnMatrix())->columnAsVector(0);
     }
 
+    /**
+     * Convolve this matrix with another matrix.
+     * 
+     * @param  \Rubix\Tensor\Matrix  $b
+     * @throws \InvalidArgumentException
+     * @return self
+     */
+    public function convolve(Matrix $b) : self
+    {
+        list($m, $n) = $b->shape();
+
+        if ($m > $this->m or $n > $this->n) {
+            throw new InvalidArgumentException('Matrix B cannot be'
+                . ' larger than Matrix A.');
+        }
+
+        $p = intdiv($m, 2);
+        $q = intdiv($n, 2);
+
+        $c = [];
+
+        for ($i = 0; $i < $this->m; $i++) {
+            $temp = [];
+
+            for ($j = 0; $j < $this->n; $j++) {
+                $sigma = 0;
+
+                foreach ($b as $k => $rowB) {
+                    foreach ($rowB as $l => $valueB) {
+                        $x = $i + ($p - (int) $k);
+                        $y = $j + ($q - (int) $l);
+
+                        if ($x < 0 or $x >= $this->n or $y < 0 or $y >= $this->m) {
+                            continue 1;
+                        }
+
+                        $sigma += $this->a[$x][$y] * $valueB;
+                    }
+                }
+
+                $temp[] = $sigma;
+            }
+
+            $c[] = $temp;
+        }
+
+        return self::quick($c);
+    }
+
     /**
      * Calculate the row echelon form (REF) of the matrix. Return the matrix
      * and the number of swaps in a tuple.

src/Vector.php

@@ -556,6 +556,35 @@ public function cross(Vector $b) : self
         return self::quick($c);
     }
 
+    /**
+     * Convolve this vector with another vector.
+     * 
+     * @param  \Rubix\Tensor\Vector  $b
+     * @throws \InvalidArgumentException
+     * @return self
+     */
+    public function convolve(Vector $b) : self
+    {
+        if ($b->size() > $this->n) {
+            throw new InvalidArgumentException('Vector B cannot be'
+                . ' longer than Vector A.');
+        }
+
+        $c = [];
+
+        for ($i = 0; $i < $this->n; $i++) {
+            $sigma = 0;
+
+            foreach ($b as $j => $valueB) {
+                $sigma += ($this->a[$i - (int) $j] ?? 0) * $valueB;
+            }
+
+            $c[] = $sigma;
+        }
+
+        return self::quick($c);
+    }
+
     /**
      * Project this vector on another vector.
      * 

tests/MatrixTest.php

@@ -564,6 +564,32 @@ public function test_matmul()
         $this->assertEquals($y, $z->asArray());
     }
 
+    public function test_convolve()
+    {
+        $input = Matrix::quick([
+            [3, 27, 66, 29, 42, 1],
+            [5, 9, 15, 42, 45, 16],
+            [91, 67, 49, 22, 66, 5],
+            [5, 1, 4, 9, 8, 6, 2],
+            [22, 16, 18, 19, 21, 25],
+            [6, 9, 69, 5, 2, 33, 35],
+        ]);
+
+        $z = $input->convolve($this->a);
+
+        $y = [
+            [254, 792, 1565, 811, 499, 195],
+            [540, 2311, 711, 2350, -766, 409],
+            [1356, 1083, 1304, 992, 478, -584],
+            [831, 164, -75, 1225, 21, -747],
+            [392, 1670, -566, 1114, 1036, -412],
+            [290, 429, 889, 69, 347, 20],
+        ];
+
+        $this->assertInstanceOf(Matrix::class, $z);
+        $this->assertEquals($y, $z->asArray());
+    }
+
     public function test_multiply_matrix()
     {
         $z = $this->a->multiply($this->c);

tests/VectorTest.php

@@ -331,6 +331,16 @@ public function test_cross()
         $this->assertEquals($y, $z->asArray());
     }
 
+    public function test_convolve()
+    {
+        $z = $this->a->convolve($this->c);
+
+        $y = [-60.0, 2.5, 259.0, -144., 40.5, 370.1, 462.20000000000005, 9.999999999999886];
+
+        $this->assertInstanceOf(Vector::class, $z);
+        $this->assertEquals($y, $z->asArray());
+    }
+
     public function test_project()
     {
         $z = $this->a->project($this->b);