Arithmetic Operations

[micuat]

Implemented + - * between matrices and * / between a matrix and a scalar. I'm new to D and not sure this is a good implementation, and also I didn't make unittests.

[kyllingstad]

Hi, and thanks for making this pull request!

However, I am sorry to say that I will not merge it. I will explain why, and I hope this rejection does not deter you from making further pull requests later. Contributions are always welcome. :-)

SciD actually had a Matrix type once, which contained, among other things, the functionality you have implemented here. I later removed it again and replaced it with the MatrixView type, which is (by design) a minimal matrix abstraction – simply a matrix-like view of a one-dimensional array, without any bells and whistles.

Here's the reason: If we add matrix operations, we should implement them properly, and that is a lot of work and requires great skill.

Firstly, when it comes to linear algebra, people have extremely high expectations with regards to performance. This is why I didn't implement the scid.linalg algorithms myself, but rather call on LAPACK functions to perform the actual calculations. That way, people can link to whichever LAPACK implementation they prefer, including one that is tuned for the hardware they use.

Secondly, in a language like D, people nowadays expect to be able to write stuff like

Matrix p, q, r;
double a, b, c;
...
auto m1 = a*p + b*q;
auto m2 = c*m1*r;
writeln(m2);

and not actually have any calculations performed until the last line, and then have several operations (perhaps all) performed at once. This significantly reduces the number of times it is necessary to loop over the elements of each matrix, and obviates the need for intermediate memory allocations.

(This is done using expression templates; see for example the C++ library Eigen. Cristi Cobzarenco and David Simcha started on something like this for D, building on SciD, a couple of years ago, but never finished it.)

If we add matrix operations, this is the standard against which we will be measured, and in my opinion it is therefore better to leave them out completely rather than have something half-baked with subpar performance. Also note that addition and subtraction, as well as element-wise multiplication and division, can be performed directly on the underlying arrays using D's array operations:

double[] p, q, r;
// Allocate and fill arrays
p[] = q[] + r[];

[micuat]

Thanks for the details. I will temporarily extend MatrixView to achieve them.

By the way, I'm thinking implementing SVD related functions into linalg. Can I do it and open a new pull request?

[kyllingstad]

Sure! You should look at the implementations of other functions in std.linalg, so that the SVD functions are written in the same style. Specifically, you should:

• Use the LAPACK SVD functions (unless there is a very good reason to implement it in D).
• Support 32- and 64-bit real and complex floating-point types.
• If possible, provide separate underscore-suffixed functions when it is necessary to modify input matrices (see std.linalg docs).
• Use scid.core.memory for temporary workspace memory management.
• Add unittests. ;-)

Are you familiar with LAPACK? Here are some links:

• LAPACK homepage
• LAPACK User's Guide: Singular Value Decomposition

[kyllingstad]

@micuat:

Thanks for the details. I will temporarily extend MatrixView to achieve them.

I will be happy to give you some tips with regard to the code you submitted here. You can consider them guidelines for further contributions. :-)

[micuat]

Sorry for posting on an old topic, but does this make sense for avoiding intermediate instances?

    void delegate(ref T[]) opBinary(string op)(MatrixView!(T) rhs)
    in
    {
        static if (op == "+" || op == "-") {
            assert (rows == rhs.rows && cols == rhs.cols);
        }
    }
    body
    {
        static if (op == "+") {
            return delegate void(ref T[] output) { output = array.dup; output[] += rhs.array[]; };
        } else static if (op == "-") {
            return delegate void(ref T[] output) { output = array.dup; output[] -= rhs.array[]; };
        } else static assert(0, "Operator "~op~" not implemented");
    }

    void opAssign(void delegate(ref T[]) rhs)
    {
        rhs(array);
    }

I know this itself is not enough for a full implementation, but I just want to know if this is conceptually right or not for implementing the operators.

[kyllingstad]

No, I'm pretty sure this will be less efficient than the naïve implementation, because:

1. You're still copying the matrix once per operation with the .dup.
2. The use of delegate could create a closure (i.e. a copy of the delegate's context; array in this case) if you're not careful, plus a delegate call is always indirect (i.e. the pointer needs to be dereferenced before the function is called).

For maximum performance, the full operations have to be constructed at compile time, using templates and mixins. See the "expression templates" Wikipedia article for what I mean. You should probably have a look at Cristi Cobzarenco's scid fork, in particular the scid.matrix module and the scid.ops package.