Vector multiplication on a 2D array using Accelerate in Swift

Question

I have the following NumPy example which does element-wise vector multiplication on a 2D array:

import numpy as np

a = np.array([[3, 5, 8, 7, 2],
              [1, 4, 9, 1, 3],
              [2, 7, 3, 1, 2],
              [5, 4, 9, 1, 6],
              [4, 2, 9, 6, 7]])

b = np.array([[2],
              [0.5],
              [9],
              [1],
              [18]])

c = a * b

print(c)

This prints the following:

[[  6.   10.   16.   14.    4. ]
 [  0.5   2.    4.5   0.5   1.5]
 [ 18.   63.   27.    9.   18. ]
 [  5.    4.    9.    1.    6. ]
 [ 72.   36.  162.  108.  126. ]]

I'm trying to do this calculation in Swift using Accelerate as shown below:

import Foundation
import Accelerate

let a: [Float] = [3, 5, 8, 7, 2,
                  1, 4, 9, 1, 3,
                  2, 7, 3, 1, 2,
                  5, 4, 9, 1, 6,
                  4, 2, 9, 6, 7]

let b: [Float] = [2,
                  0.5,
                  9,
                  1,
                  18]

var c: [Float] = Array(repeating: 0, count: 25)
vDSP_vmul(a, 1, b, 1, &c, 1, 25)

print(c)

The Swift code prints the following which is not correct:

[6.0, 2.5, 72.0, 7.0, 36.0, 0.0, 0.0, 0.0, 1.104904e-39, 4.7331654e-30, 3.7740115e-24, 6.737052e-37, 1.806498e-35, 1.19422e-39, 5.1722454e-26, 5.9787867e-25, 2.7286e-41, 0.0, 0.0, 0.0, 1.6543612e-24, 0.0, 0.0, 5.533857e-39, 3.73e-43]

Accelerate works with flatten arrays whereas NumPy supports 2D arrays. So I'm not sure how to multiply the 2D array with a vector in Accelerate. Is there a different Accelerate function that I should use for this type of multiplication?

Answer 1

Just to throw it out there, here’s an alternative: Given that there is no built-in vector-element-wise product function, rather than using a vector, you could use a diagonal matrix and then do a simple matrix product (e.g., cblas_sgemm):

import Accelerate

// render vector as a diagonal matrix

var a: [Float] = [2,   0,   0,   0,   0,
                  0,   0.5, 0,   0,   0,
                  0,   0,   9,   0,   0,
                  0,   0,   0,   1,   0,
                  0,   0,   0,   0,  18]

let rowsA: Int32 = 5
let columnsA: Int32 = 5

var b: [Float] = [3,   5,   8,   7,   2,
                  1,   4,   9,   1,   3,
                  2,   7,   3,   1,   2,
                  5,   4,   9,   1,   6,
                  4,   2,   9,   6,   7]

let rowsB: Int32 = 5
let columnsB: Int32 = 5

assert(columnsA == rowsB)

var c: [Float] = Array(repeating: 0, count: Int(rowsA * columnsB)) // Result matrix

// Perform matrix multiplication using cblas_sgemm

cblas_sgemm(
    CblasRowMajor, // Matrix layout
    CblasNoTrans,  // Transpose option for matrix A
    CblasNoTrans,  // Transpose option for matrix B
    rowsA,         // Number of rows in matrices A and C
    columnsB,      // Number of columns in matrices B and C
    columnsA,      // Number of columns in matrix A; number of rows in matrix B
    1,             // Scaling factor for the product of matrices A and B
    &a,            // Pointer to matrix A
    columnsA,      // Leading dimension of matrix A
    &b,            // Pointer to matrix B
    columnsB,      // Leading dimension of matrix B
    0,             // Scaling factor for matrix C
    &c,            // Pointer to the result matrix
    columnsB       // Leading dimension of the result matrix
)

// Print the result matrix

for row in 0..<rowsA {
    for column in 0..<columnsB {
        print(c[Int(row * columnsB + column)], terminator: " ")
    }
    print()
}

That yields:

  6.0  10.0  16.0  14.0   4.0 
  0.5   2.0   4.5   0.5   1.5 
 18.0  63.0  27.0   9.0  18.0 
  5.0   4.0   9.0   1.0   6.0 
 72.0  36.0 162.0 108.0 126.0 

Answer 2

To multiply this as a column, you need to use a stride. If you're taking the trouble to use Accelerate, then you probably also want to avoid extra allocations. You can use Array(unsafeUninitializedCapacity:) to fill in all the data like this:

let c = Array(unsafeUninitializedCapacity: a.count) { output, initializedCount in
    a.withUnsafeBufferPointer { abuffer in
        for i in b.indices {
            vDSP_vmul(abuffer.baseAddress! + i, b.count, // stride through a by the column size
                      b, 1,                              // stride one-by-one through b
                      output.baseAddress! + i, b.count,  // stride through the output by the column size
                      vDSP_Length(b.count))
        }
        initializedCount = abuffer.count
    }
}

// [6.0, 10.0, 16.0, 14.0, 4.0, 0.5, 2.0, 4.5, 0.5, 1.5, 18.0, 63.0, 27.0, 9.0, 18.0, 5.0, 4.0, 9.0, 1.0, 6.0, 72.0, 36.0, 162.0, 108.0, 126.0]