Broadcasting rules

[Feras Salim]

Wondering what are the broadcasting rules as the following seems to work fine in xtensor numpy and matlab

    const double K =100;
    const double r = 0.02;
    const double div = 0.05;
    const double sigma_sq = 0.2*0.2;
    const double L1 = -5;
    const double L2 = 1.5;
    const unsigned long N = ceil((L2-L1)/dx+1);
    const unsigned long M = 1095;

    blaze::DynamicMatrix<double,blaze::rowMajor> x = expand(linspace<rowVector>( N, L1, L2),1);
    blaze::DynamicMatrix<double,blaze::columnMajor> tau = expand(linspace<columnVector>( M, 0., T),1);
    cout << x.rows() << " X " << x.columns() << endl;
    cout << tau.rows() << " X " << tau.columns() << endl;
    blaze::DynamicMatrix<double> S = (K*exp(x-(r-div-0.5*sigma_sq)*tau));

1 X 1076
1095 X 1
unknown file: Failure
C++ exception with description "Matrix sizes do not match" thrown in the test body.

The expected result is 1095x1076 matrix.

A couple of side questions, if I do it without expand it also doesn’t work which is in line with most other libraries:

1. Do you treat vectors as single row/column matrices when broadcasting or not (the numpy way is to not)
2. Do you expose any debugging facilities to print out a nicer version of dimensions and truncated matrix when too large (https://xtensor.readthedocs.io/en/latest/api/xio.html)
3. Is there any meaningful overhead in calling expand with runtime params. For instance with C++17 I can’t seem to do blaze::DynamicMatrix<double,blaze::rowMajor> x = expand<1>(linspace<rowVector>( N, L1, L2));
4. Is there a better means to ask some of these higher level questions like gitter or mailing list or the preferred way is tickets?

Thanks!

‌

[Klaus Iglberger]

Hi Feras!

1. No, also in Blaze we do not. The problem in your code, however, is not the broadcasting rules, but the precedence of operations. Essentially you are trying to subtract a 1095x1 matrix from a 1x1076 matrix (in short x-factor*tau, multiplication before subtraction).

2. Since the user requirements on “nice” are very diverse, you’ll have to provide your own pretty print functionality. We are open for suggestions, though.

3. There is a small overhead, which in your case is negligible, though, because your “runtime” parameter is also known at compile time. Your C++17 example compiles and works fine with GCC-9, Clang-9.0 and MSVC 2017. In case you encounter a real compilation issue, please create an issue with the exact information under which conditions the code does not compile.

4. We are currently setting up a new forum for this kind of question (see issue #170). Hence at this point the only two ways are to create issues or to use our old forum (https://groups.google.com/forum/#!forum/blaze-forum).

Best regards,

Klaus!

[Feras Salim]

Hi Klaus,

Thanks for the quick reply.

Re 1. I’m basically trying to reproduce this : https://repl.it/@feribg/broadcasting

Re 2. I guess something as simple as not printing the entire matrix when it’s too large and truncating with … row and column wise with a configurable truncation threshold would suffice.

Re 3. It works fine (just triggers a warning with Clang 11 on OSX) for the first case but the second case fails, I assumed row and column vectors are orthogonal so when you call expand the rows and cols should be mirror of eachother, correct?

blaze::DynamicMatrix<double,blaze::columnMajor> x = expand<1>(linspace<columnVector>( 10, 1, 10));
EXPECT_EQ(10, x.rows());
EXPECT_EQ(1, x.columns());

#Fails, dims are the same as above
blaze::DynamicMatrix<double,blaze::rowMajor> y = expand<1>(linspace<rowVector>( 10, 1, 10));
EXPECT_EQ(1, y.rows());
EXPECT_EQ(10, y.columns());

warning: use of function template name with no prior declaration in function call with explicit template arguments is a C++2a extension [-Wc++2a-extensions]
    blaze::DynamicMatrix<double,blaze::rowMajor> x = expand<1>(linspace<rowVector>( 1, 1, 10));

‌

Re 4. Awesome, let me know what you think is the best way, here or the google group. I see the google group is quite stale, so maybe just adding a question label in the issue tracker would be fine?

[Klaus Iglberger]

Hi Feras!

1. I understand your question. In Blaze you cannot subtract a column vector from a row vector or vice versa. Both vectors need to be the same kind of vector. If you need to convert between row and column vectors, you can use the trans() function.
2. Thanks!
3. When I copy-and-paste the two lines for x and y, I get the following output:

blaze::DynamicMatrix<double,blaze::columnMajor> x = expand<1>(linspace<columnVector>( 10, 1, 10));
std::cerr << "x: " << rows(x) << "x" << columns(x) << "\n";
std::cerr << x << "\n\n";

blaze::DynamicMatrix<double,blaze::rowMajor> y = expand<1>(linspace<rowVector>( 10, 1, 10));
std::cerr << "y: " << rows(y) << "x" << columns(y) << "\n";
std::cerr << y << "\n\n";

x: 10x1
(  1 )
(  2 )
(  3 )
(  4 )
(  5 )
(  6 )
(  7 )
(  8 )
(  9 )
( 10 )

y: 1x10
( 1 2 3 4 5 6 7 8 9 10 )

It sounds like this is exactly what you would expect. x is a 10x1 matrix, y is a 1x10 matrix.

4. Both is fine from my point of view. The old forum isn’t used anymore (hence issue #170), but should still serve well for asking questions.

Best regards,

Klaus!

[Feras Salim]

Thanks Klaus!

Re 1. Is that the case for all arithmetic operations (and maybe not just arithmetic) or just addition and subtraction. What would be the equivalent Blaze code to achieve that but still preserving contiguous data? I assume trans is just a view over the same row major data that’s being used as a column. In the current case x and tau are no longer vectors but matrices, so it seems a common use case.

Best,

Feras

[Klaus Iglberger]

Hi Feras!

The wiki describes the rules for the different arithmetic operations. Here is a short summary for the most common operations:

blaze::DynamicVector<int,columnVector> c1, c2;
blaze::DynamicVector<int,rowVector>    r1, r2;

blaze::DynamicMatrix<int,rowMajor> R1, R2;
blaze::DynamicMatrix<int,columnMajor> C1, C2;

// Vector Addition
c1 + c2;  // Vector addition of two column vectors (vector size must match)
r1 + r2;  // Vector addition of two row vectors (vector size must match)
c1 + r2;  // Outer sum of a column vector and a row vector
c1 + 99;  // Vector/scalar addition (both for column and row vectors)
99 + c2;  // Scalar/vector addition (both for column and row vectors)

// Vector Subtraction
c1 - c2;  // Vector subtraction of two column vectors (vector size must match)
r1 - r2;  // Vector subtraction of two row vectors (vector size must match)
c1 - r2;  // Outer difference of a column vector and a row vector
c1 - 99;  // Vector/scalar subtraction (both for column and row vectors)
99 - c2;  // Scalar/vector subtraction (both for column and row vectors)

// Vector multiplication
c1 * c2;  // Componentwise vector multiplication of two column vectors (vector size must match)
r1 * r2;  // Componentwise vector multiplication of two row vectors (vector size must match)
c1 * r2;  // Outer product of a column vector and a row vector
r1 * c2;  // Inner product of a row vector and a column vector (vector size must match)
c1 * 99;  // Scalar multiplication (both for column and row vectors)
99 * c2;  // Scalar multiplication (both for column and row vectors)

// Matrix/vector multiplication
R1 * c1;  // Matrix/vector multiplication (row-major matrix and column vector; number of columns and vector size must match)
C1 * c1;  // Matrix/vector multiplication (column-major matrix and column vector; number of columns and vector size must match)

// Vector/matrix multiplication
r1 * R1;  // Vector/matrix multiplication (row vector and row-major matrix; vector size and number of rows must match)
r1 * C1;  // Vector/matrix multiplication (row vector and column-major matrix; vector size and number of rows must match)

// Matrix addition
R1 + R2;  // Matrix/matrix addition (two row-major matrices; number of rows and columns must match)
R1 + C2;  // Matrix/matrix addition (row-major and column-major matrix; number of rows and columns must match)
C1 + R2;  // Matrix/matrix addition (column-major and row-major matrix; number of rows and columns must match)
C1 + C2;  // Matrix/matrix addition (two column-major matrices; number of rows and columns must match)
R1 + 99;  // Matrix/scalar addition (both for row-major and column-major matrices)
99 + R2;  // Matrix/scalar addition (both for row-major and column-major matrices)

// Matrix subtraction
R1 - R2;  // Matrix/matrix subtraction (two row-major matrices; number of rows and columns must match)
R1 - C2;  // Matrix/matrix subtraction (row-major and column-major matrix; number of rows and columns must match)
C1 - R2;  // Matrix/matrix subtraction (column-major and row-major matrix; number of rows and columns must match)
C1 - C2;  // Matrix/matrix subtraction (two column-major matrices; number of rows and columns must match)
R1 - 99;  // Matrix/scalar subtraction (both for row-major and column-major matrices)
99 - R2;  // Matrix/scalar subtraction (both for row-major and column-major matrices)

// Matrix multiplication
R1 * R2;  // Matrix/matrix product (two row-major matrices; number of rows of R1 must match number of columns of R2)
R1 * C2;  // Matrix/matrix product (row-major and column-major matrix; number of rows of R1 must match number of columns of C2)
C1 * R2;  // Matrix/matrix product (column-major and row-major matrix; number of rows of C1 must match number of columns of R2)
C1 * C2;  // Matrix/matrix product (two column-major matrices; number of rows of C1 must match number of columns of C2)
R1 * 99;  // Scalar multiplication (both for row-major and column-major matrices)
99 * R2;  // Scalar multiplication (both for row-major and column-major matrices)

In your example you are trying to subtract two matrices with different dimensions (x-tau). This will cause a runtime failure. You’ll be able to subtract them if the dimensions match (e.g. x-trans(tau)). There is no implicit broadcast, however. For that purpose you can use the outer difference between two vectors.

I hope this helps,

Best regards,

Klaus!

[Feras Salim]

Thanks Klaus,

The outer product works fine for vector - vector (just need to change the order in the operation because it seems to matter), but I still can’t find a way to do the numpy way of doing broadcasting, IE when there’s a unit axis, expand to the size of the matrix. Seems I have to do it explicitly in Blaze, is that the case and the recommended way of doing it, the expand operator doesn’t allocate, so the overhead should be minimal if used like that ?

//S shape = 1095x1076

//works
auto M = 1095;
auto tau = linspace<columnVector>( M, 0., T);
bound = (S-K) % exp(r*(expand(tau,S.columns())));

//doesn't work
bound = (S-K) % exp(r*tau)

//doesn't work
bound = (S-K) % exp(r*(expand(tau,1)));

Best,
Feras

[Klaus Iglberger]

Hi Feras!

You are correct, there is no implicit broadcast operation in Blaze. All operations in Blaze perform a specific action and don’t perform different actions based on the state of the input. Also, the expand() operation never allocates, but returns a proxy object (a so-called expression object) that pretends to be the accordingly sized matrix.

Best regards,

Klaus!