1

I am learning about template meta programming and expression templates in C++ right now, so as an exercise, I am creating a linear algebra library to practice the concepts I am learning.

So far, my library has a complete list of non-member operator overloads for all of the binary operators that can be overloaded, and has a fairly-slick interface that's easily expandable. One problem I've run into, however, is that matrix operations often have multiple variations. For example, for multiplication, there is the general matrix multiplication, the dot product, the kroenecker product, the hadamard product, the cross product, and more.

One slick way around this that is employed in Matlab is the .* operator used for hadamard multiplication (and .^, ./, etc). In this case, the Matlab language employs the . operator as a modifier for the * operator. However, I'm unaware of any mechanisms in the c++ language that allow operators to be modified like this. Are there any clean workarounds for this behavior?

Here are some things I've thought about already:

  • operator overloads allow extra template parameters. However, I'm not entirely sure how to take advantage of this in this context. For instance, something that might be nice (though, in practice, I am not sure there is a valid syntax to achieve this):
template<typename lhs_t, typename rhs_t, typename op_t = Gemm>
auto operator*(lhs_t lhs, rhs_t rhs)
{
    ...
}

// Operator template specializations ...

int main()
{
    Matrix<double, 2, 2> mat1(1.0, 2.0, 3.0, 4.0);
    Matrix<double, 2, 2> mat2(1.0, 2.0, 3.0, 4.0);

    mat1 * mat2; // Works
    mat1 *<Hadamard> mat2; // Error!  Syntax????
}
  • Using SFINAE/Concepts/if constexpr and traits to modify binary expression types or wrap binary expression types. Syntax:
Hadamard(mat1 * mat2); // Hadamard wraps or modifies binary expression created by operator*
                       // SFINAE or Concepts used to select the correct routine based on the trait set
  • Create a free binary function. Possible syntaxes:
Hadamard<Multiplication>(mat1, mat2);
Hadamard_Multiplication(mat1, mat2);
  • Using member functions. Syntax:
mat1.hadamard_multiplication(mat2);

None of these seem to have syntax quite as elegant as Matlab's:

mat1 .* mat2;

Are there any techniques available to come close to an operator modifier syntax that I can consider? Or any general techiques to make the usage syntax less verbose? If not, is there any notion that something could be included in future versions of C++ that may be of any use here?

  • 2
    `.*` is one of the ways to access a member via a [pointer to member](https://en.cppreference.com/w/cpp/language/operator_member_access). – 1201ProgramAlarm Jan 26 '20 at 22:06
  • 1
    You can only modify an existing operator, or use a function. `hadamard(mat1, mat2);` would be a good style . – M.M Jan 26 '20 at 22:06
  • 1
    You cannot create new operators in C++. – Sam Varshavchik Jan 26 '20 at 22:17
  • I understand why `.*` and `.` cannot be overloaded, but it is a shame. Really, any form of compounded operators would be useful. – Christopher Mauer Jan 26 '20 at 22:31
  • 1
    You could overload `operator->*`, or you could overload `operator()` to return some wrapper object, which itself overloads `operator*(Mat const&)` (allowing either `mat1 ->* mat2` or `mat1() * mat2`), though I'd argue that often verbosity is not bad. E.g., [Eigen](http://eigen.tuxfamily.org/) allows these syntaxes: `mat1.array() * mat2.array()` or `mat1.cwiseProduct(mat2)` – chtz Jan 28 '20 at 13:26
  • 1
    "None of these seem to have syntax quite as elegant as Matlab's", but Hadamard(mat1 * mat2) appears to be much more readable. – StefanKssmr Jan 29 '20 at 06:54
  • I think that depends if you are a software developer putting on the hat of a mathematician, or a mathematician putting on the hat of a software developer. I mostly work with the latter, who would disagree with you. But I do personally agree with you. – Christopher Mauer Jan 29 '20 at 12:44
  • 1
    With the `Hadamard(mat1 * mat2)` syntax, how would you express, e.g., `(A*B) .* C`? You could write `Hadamard(A*B*C)` and for `A.*B.*C` write `Hadamard(Hadamard(A*B) * C)`, but I would not consider that "less verbose" than any alternative. – chtz Jan 29 '20 at 13:58
  • I like the `hadamard(mat1*mat2)` syntax because it negates needing to specify that it's a multiplication (as opposed to hadamard exponentiation/division/etc), but your example of `hadamard(A*B*C)` is interesting. I think it would expand to: hadamard(multiplication(multiplication(A,B),C)), which is not what is desired. Furthermore, this would be a completely legal under that syntax, and may be a bit confusing. I had not thought of that. – Christopher Mauer Jan 29 '20 at 14:08

1 Answers1

0

Here is the realization that I came to with respect to this:

  1. Multiple syntaxes can be supported. I do not need to pick one, I just need to implement one way and add abstractions over the top mimick the other syntaxes if so desired.
  2. Hadamard multiplication (and division, etc) are the default multiplication style for arrays, while GEMM is the default multiplication mode for matrices.

With those two items, I chose to implement an n-dimensional array class (since CRTP base classes and free functions are used for the implementation, this was simply a matter of providing an empty struct with the necessary traits/alias declarations).

The following situations can then result in hadamard multiplications:

  1. Matrix A has the same dimensions as Matrix B, and neither matrix is a square matrix, it is then clear that gemm is not valid, but hadamard multiplication is. Therefore, it makes sense to use a concept to override this behavior. in other words:

    // Results in hadamard multiplication
    template<object_type lhs_t, object_type rhs_t> requires 
    (is_same_dimensions_v<lhs_t, rhs_t> && !is_special_mult_supported_v<lhs_t, rhs_t>)
    constexpr auto operator*(lhs_t&& lhs, rhs_t&& rhs) noexcept -> mul_expr<lhs_t, rhs_t>
    {
        return { std::forward<lhs_t>(lhs), std::forward<rhs_t>(rhs) };
    }
    
    // Results in general matrix multiplication
    template<typename lhs_t, typename rhs_t> requires is_gemm_supported_v<lhs_t, rhs_t>
    constexpr auto operator*(lhs_t&& lhs, rhs_t&& rhs) noexcept -> gemm_expr<lhs_t, rhs_t>
    {
        return { std::forward<lhs_t>(lhs), std::forward<rhs_t>(rhs) };
    }
    
  2. If a gemm expression is assigned to an array, and a hadamard multiplication is valid, it is implicitly converted to a hadamard multiplication.

    // implicit conversion to hadamard multiplication
    array = a * b;
    
  3. Gemm expressions can be cast to hadamard expressions explicitly.

    // explicitly create a hadamard multiplication expression
    auto c = static_cast<mul_expr>(a * b);
    
  4. Gemm expressions can be cast to arrays explicitly, resulting in a hadamard multiplication

    // explicitly create an array using hadamard multiplication
    auto c = static_cast<array>(a * b);
    
  5. In the case of mixed array/matrix types where both hadamard and gemm are supported, the left-hand side type is selected to prevent ambiguity.

     // if a is an array, and b is a matrix, then c is a mult_expr
     // if a is a matrix, and b is an array, then c is a gemm_expr
     auto c = a * b;
    

With these rules in place, then api-level abstractions can be added for clarity:

    // c = hadamard(a * b);
    template<expr_type expr_t> requires std::is_convertible_v<std::decay_t<expr_t>, mul_expr>
    constexpr auto hadamard(expr_t&& expr) noexcept
    {
        return static_cast<mul_expr>(expr);
    }

    // support c = hadamard_mult(a, b); syntax
    template<object_type lhs_t, object_type rhs_t> requires is_same_dimensions_v<lhs_t, rhs_t>
    constexpr auto hadamard_mult(lhs_t&& lhs, rhs_t&& rhs) noexcept
    {
        return hadamard(lhs * rhs);
    }

Note that I've omitted static_casts and paraphrased some of the syntax to get the idea across. The big realization to take away here is that c++ can utilize the type system to simplify syntax fairly drastically, which is where matlab differs. There are many cases where

    c = a * b;

can (and should) result in hadamard multiplication. Furthermore, simple manipulation of the type system can quite easily result in situations where function syntax is easily supported.

A big thank you to those of you in the comments above for helping me brainstorm and arrive at this conclusion. I am quite satisfied with my library as a consequence of the discussion here.