Questions around element-wise treatments & matrix multiplication syntax

[Phil]

• Does blaze support element-wise operations for matrices? The examples below unfortunately do not compile:

StaticMatrix<int, 3, 2> M{ {3, 8},{7, 6},{2, 4} };
StaticMatrix<int, 3, 2> N{ {4, 7},{5, 9},{6, 7} };

auto result_3x2_div_3x2 = M / N;  // how to make them be applied element-wise?
auto result_3x2_mul_3x2 = M * N;  // here i do NOT want matrix multiplication

‌

• Sometimes operator* is element-wise and sometimes it becomes a matrix-product. Why not reserve operator* for element-wise operations only (just like operator+ and operator-)? For consistency. And then use e.g. mult(M, v) to explicitly multiply two objects.

DynamicMatrix<double> M{ {3, 8},{8, 7},{5, 4} };
DynamicVector<double> v{ 4, 9 };
auto result1 = M * v;  // matrix product: M.v = {84, 95, 56}
auto result2 = v * v;  // elmt-wise multiplication: v[0]*v[0], v[1]*v[1]} = {16, 81}
                       // --> inconsistent use of operator*

‌

• Is it possible to add (& subtract, multiply, divide) N×M matrices and K-sized vectors, for when K == M && K != N or K == N && K != M? Obviously, the case K == M == N is ill-defined.

StaticMatrix<int, 3, 2> M{ {3, 8},{7, 6},{2, 4} };
StaticVector<int, 3> v{ 4, 6, 9 };
auto result_3x2_add_3x1 = M + v;  // <-- this line fails

‌

• How come one cannot convert a matrix-vector product back to a matrix? The result is a 3-sized vector, or in other words a 3×1 matrix.

DynamicMatrix<double> M{ {3, 8},{8, 7},{5, 4} };
DynamicVector<double> v{ 4, 9 };

auto result1                  = M * v;  // <-- works, it's just the expression
DynamicVector<double> result2 = M * v;  // <-- works, convert expr back to vec
DynamicMatrix<double> result3 = M * v;  // <-- this line fails

‌

• Is there an expression to reshape objects? Say we have a vector of size K and want to reshape it into an N×M matrix (with N*M == K of course). Same principle for reshaping N×M matrices to either K×L matrices (with N*M == K*L) or to vectors of size K (with K == N*M). For example:

DynamicVector<double> v{ 1,2,3, 7,8,9 };
auto matrix_2x3 = reshape(v, 2, 3);  // should give {{1,2,3},{7,8,9}}

‌

Happy to discuss further, and hope my questions make sense.

[Klaus Iglberger]

Hi Phil!

Does blaze support element-wise operations for matrices? The examples below unfortunately do not compile:

Yes, Blaze supports element-wise operations for matrices. In fact, most of the available matrix operations are element-wise operations. You are even able to define your own element-wise matrix operations via the map() function. With map() you are for instance able to define matrix/matrix division:

template< typename MT1
        , bool SO1
        , typename MT2
        , bool SO2 >
decltype(auto) operator/( const blaze::DenseMatrix<MT1,SO1>& lhs, const blaze::DenseMatrix<MT2,SO2>& rhs )
{
   return blaze::map( *lhs, *rhs, Div{} );
}

The Schur product is expressed via operator%:

blaze::DynamicMatrix<double>   M1( 28UL, 35UL );
blaze::CompressedMatrix<float> M2( 28UL, 35UL );

// ... Initialization of the matrices

DynamicMatrix<double> M3 = M1 % M2;

Sometimes operator* is element-wise and sometimes it becomes a matrix-product. Why not reserve operator* for element-wise operations only (just like operator+ and operator-)? For consistency. And then use e.g. mult(M, v) to explicitly multiply two objects.

A consistent interface is an important argument, but more important is an intuitive and non-surprising interface. For matrices, operator* is interpreted as the matrix product by most people and therefore chosen for this operation. Also, since the Schur product is not used in many applications (it for instance is not even considered in linear algebra proposal for the C++ standard), it uses the “special” operator%.

Most people would consider mult() the same as operator*() (i.e. _mult_iplication). Therefore using mult() for a different operation would again be surprising and non-intuitive. Since the learning effort would be the same as for the separation into operator*() and operator%(), we decided to stick to the more intuitive interface.

Is it possible to add (& subtract, multiply, divide) N×M matrices and K-sized vectors, for when K == M && K != N or K == N && K != M? Obviously, the case K == M == N is ill-defined.

No. Matrices and vectors are strictly separated in order to provide as much compile time optimization as possible. You are, however, able to extract vectors from matrices via the row() and column() functions, and you are able to expand() a vector into a matrix to express your intention. For your convenience, you are for instance able to provide a custom toMatrix() function:

template< typename VT, bool TF >
decltype(auto) toMatrix( const DenseVector<VT,TF>& vec )
{
   return expand<1UL>( vec );
}

blaze::DynamicMatrix<int> A{ { 2, 4 }, { 6, 8 } };
blaze::DynamicVector<int> v{ 1, 2 };
blaze::DynamicMatrix<int> C;

C = toMatrix( A * v );

How come one cannot convert a matrix-vector product back to a matrix? The result is a 3-sized vector, or in other words a 3×1 matrix.

As stated above, vectors and matrices are strictly separated in order to provide as much compile time optimization as possible. You can assign to a row() or column() of a matrix or expand() the result.

Is there an expression to reshape objects? Say we have a vector of size K and want to reshape it into an N×M matrix (with N*M == K of course). Same principle for reshaping N×M matrices to either K×L matrices (with N*M == K*L) or to vectors of size K (with K == N*M). For example:

A reshape() operation is a special case for matrices with NxM contiguous elements. In order to provide maximum performance, especially for tiny and small vectors and matrices, Blaze does not guarantee contiguous elements, but uses padding (see for instance the StaticMatrix documentation). Also, Blaze allows you to define submatrices, which cannot guarantee contiguous elements. As a result, an efficient reshape() operation would only work on a subclass of matrices. It would be possible to provide a general reshape() that works on all matrices, but it would be horribly slow as it would use divisions and modulo operation for every element access. Since Blaze is focused on performance, we avoid this kind of performance traps in the interface.

I hope this answers your questions,

Best regards,

Klaus!

[Phil]

Hi Klaus,

Thank you for the fast response.

As for more important is an intuitive and non-surprising interface. I realise now that my point was probably a bit subjective, because “intuitive” is usually what one knows (and there may not be a correct way anyway). I say this, because I often use numpy’s built-in algebra package (which reserves + - / * for element-wise exclusively and uses free function matmul for explicit matrix multiplication). So to me blaze’s interface is surprising. But I suppose, if one is familiar with eigen, then intuition is the other way round since eigen uses * for matrix multiplication. Maybe a bit of a long-shot, but would it be worth including a mult function in blaze, which simply mirrors operator*? That would probably benefit everyone. As you say, most people would consider mult() the same as operator*().

As for the toMatrix function. Thank you, this is very useful.

Another few questions, if you don’t mind:

• I see that IdentityMatrix is always defined as a N×N square matrix. In my work, I often deal with N×M non-square identity matrices, e.g. 2×3 would be {{1,0,0}, {0,1,0}}. Would it be possible to add support for this? Happy to open a separate issue (if that would help track it).
• Does blaze support swapping the elements of specific columns/rows? I tried blaze::swap(column(M, i), column(M, j)) to swap column i and j but it does not compile. In fact, would std::swap(column(M, i), column(M, j)) do the desired job?
• Is it possible to transform a Vector to a DiagonalMatrix? Say we have DynamicVector<double> v = {2,3,4}, can we do DiagonalMatrix<double> D(v) which then represents{{2,0,0},{0,3,0},{0,0,4}}? It doesn’t have to be an explicit matrix type, an expression that represents the vector’s internal data as a matrix would be enough.
• Blaze provides maxNorm (or equivalently linfNorm) as max(abs(*dm)). Would it be possible to add minNorm as well, which would be de defined as min(abs(*dm))? See e.g. Infinity and Negative Infinity Norm of a Vector.
• Functions pow(M, k) and pow(v, k) return M^k and v^k element-wise (for both vectors & matrices). But pow2(M) does a matrix multiplication M.M (and thus fails for non-square matrices), whereas pow(v, k)does an element-wise calculation. So pow(M,2) ≠ pow2(M) but pow(v,2) = pow2(v). Is this intended? I guess the problem is related to pow2 being defined as a * a (which again, differs in meaning for vectors & matrices).

Thanks again for your helpful response.

[Klaus Iglberger]

Hi Phil!

“intuitive” is usually what one knows

100% agreement.

As for the toMatrix function. Thank you, this is very useful.

This may indeed be something for the next Blaze release.

I see that IdentityMatrix is always defined as a N×N square matrix. In my work, I often deal with N×M non-square identity matrices, e.g. 2×3 would be {{1,0,0}, {0,1,0}}. Would it be possible to add support for this? Happy to open a separate issue (if that would help track it).

The Blaze definition of an identity matrix is based on the Wikipedia definition. We have thought about non-square identity matrices, but found that they would not be unique (e.g. a 2x3 identity matrix could be {{1,0,0},{0,1,0}} or {{0,1,0},{0,0,1}}.

I suggest that you implement a custom function for your purposes:

template< typename T >
decltype(auto) id( size_t m, size_t n )  // Other name: eye
{
   const size_t N( max(m,n) );
   return submatrix( blaze::IdentityMatrix<T>{N}, 0UL, 0UL, m, n );
}

// Creates the non-square identity matrix
//   ( 1 0 0 )
//   ( 0 1 0 )
const auto I = id<int>( 2UL, 3UL );

Does blaze support swapping the elements of specific columns/rows? I tried blaze::swap(column(M, i), column(M, j)) to swap column i and j but it does not compile. In fact, would std::swap(column(M, i), column(M, j)) do the desired job?

As the wiki explains, this is currently not possible via the swap() function. You can swap via the std::swap_ranges() algorithm, though:

blaze::DynamicMatrix<int> A
      { { 1, 2, 3 }
      , { 4, 5, 6 }
      , { 7, 8, 9 } };

auto c0 = column(A,0);
auto c2 = column(A,2);

std::swap_ranges( c0.begin(), c0.end(), c2.begin() );

Admittedly this feels a little clumsy and inconvenient. I didn’t realize this till now and will see what we can do about it.

Is it possible to transform a Vector to a DiagonalMatrix? Say we have DynamicVector<double> v = {2,3,4}, can we do  DiagonalMatrix<double> D(v) which then represents{{2,0,0},{0,3,0},{0,0,4}}? It doesn’t have to be an explicit matrix type, an expression that represents the vector’s internal data as a matrix would be enough.

Unfortunately not. The best alternative that Blaze currently provides is this:

blaze::DiagonalMatrix<DynamicMatrix<int>> D( 3UL, 3UL );
blaze::DynamicVector<int> v{ 2, 3, 4 };

// Set the diagonal of D to v
diagonal(D) = v;

Blaze provides maxNorm (or equivalently linfNorm) as max(abs(*dm)). Would it be possible to add minNorm as well, which would be de defined as min(abs(*dm))? See e.g. Infinity and Negative Infinity Norm of a Vector.

This would be a nice addition for the next release. Until then, you can define the norm yourself:

template< typename MT  // Type of the dense matrix
        , bool SO >    // Storage order
decltype(auto) minNorm( const DenseMatrix<MT,SO>& dm )
{
   return min( abs( *dm ) );
}

Functions pow(M, k) and pow(v, k) return M^k and v^k element-wise (for both vectors & matrices). But pow2(M) does a matrix multiplication M.M (and thus fails for non-square matrices), whereas pow(v, k)does an element-wise calculation. So pow(M,2) ≠ pow2(M) but  pow(v,2) = pow2(v). Is this intended? I guess the problem is related to pow2 being defined as a * a (which again, differs in meaning for vectors & matrices).

I wasn’t aware of this inconsistency till now. Thanks for pointing this out. Let me check if changing the behavior of pow2() for matrices would have any ill effect. If not, we could adapt this to the expected element-wise behavior.

‌

Thanks for pointing out so many important points. We’ll look into swapping rows/columns, setting the diagonal of a diagonal matrix via a vector, add a minNorm() function, and check the behavior of pow2(). That should keep us busy on the weekend

Best regards,

Klaus!

[Phil]

100% agreement.

Do you think it may be worth including a mult function, which is synonymous to operator* (for matrices)? I could write my own, but if possible, I prefer as little customisation as possible (or else it becomes a possible maintenance issue down the line).

We have thought about non-square identity matrices, but found that they would not be unique. (e.g. a 2x3 identity matrix could be {{1,0,0},{0,1,0}} or {{0,1,0},{0,0,1}}

I have personally never seen the second variant used. But you are right, they do seem to exist. How about having a interface similar to https://numpy.org/doc/stable/reference/generated/numpy.eye.html? The signatures could be:

IdentityMatrix(size_t n);                     // Square n×n matrix, e.g.     {{1,0,0},{0,1,0},{0,0,1}} for n=3
IdentityMatrix(size_t n, size_t m);           // Non-Square n×m matrix, e.g. {{1,0,0},{0,1,0}}}        for n=2, m=3 (effectively defaulting to k=0)
IdentityMatrix(size_t n, size_t m, size k);   // Non-Square n×m matrix, e.g. {{0,1,0},{0,0,1}}}        for n=2, m=3, k=1

‌

We’ll look into swapping rows/columns, setting the diagonal of a diagonal matrix via a vector, add a minNorm() function, and check the behavior of pow2()

Thanks a lot for looking at all these points!

[Ecolss Logan]

Thanks for asking the above questions, Phil!

I just had some of the same questions, and thank Klaus for the answers.
I was looking through wiki and resolved issues, trying to find what the native support for matrix/matrix division, just found out the map solution.

[Klaus Iglberger]

Hi Phil!

With the next push you’ll get the alternate interface functions add(), sub(), mult(), div(), and schur() and the minNorm() function. The other features are still in progress.

Best regards,

Klaus!

[Phil]

Thanks a lot, Klaus.

When you say “alternate interface functions add(), sub(), mult(), div()”, are you referring to element-wise operations (for objects with matching shapes), with mult and schur being the same? Or mult being the same as operator* (matrix multiplication)?

‌

[Klaus Iglberger]

add() will be the alternative to operator+(), sub() the alternative to operator-(), mult() the alternative for operator*(), div() the alternative to operator/() and schur() the alternative to operator%().

[Phil]

Perfect, thank you.