Date: Mon, 1 May 2023 17:43:26 -0400

> On May 1, 2023, at 5:37 PM, Mark Hoemmen via Lib-Ext <lib-ext_at_[hidden]> wrote:

>

> Greetings and thanks for your hard work!

>

> I've made it a point to limit my comments about P1385. It's a

> different design that serves a different set of users. That being

> said, I'd like to comment on outer products and rank-1 updates.

>

> On Mon, May 1, 2023 at 5:55 AM Guy Davidson via SG19

> <sg19_at_[hidden] <mailto:sg19_at_[hidden]>> wrote:

>>

>> Hello everyone

>>

>> I'm just putting together a first pass at the wording for P1385, A proposal to add linear algebra support to the C++ standard library. If you look at the latest revision you will infer that at Kona in November and at Issaquah in February I addressed SG6 and LEWG about withdrawing the vector class entirely and simply offering a matrix class, where a vector is a special case of a matrix, with a single row or column. There were no objections to this approach.

>>

>> While there were no objections raised in the meeting, others have come in, and I want to use the reflectors to gather opinion about the matter. The heart of the problem is: what does the vector product signify? Is it an inner or outer product? Is vector orientation significant?

>>

>> With my mathematician's hat on, multiplying a row vector by a column vector is an inner product, yielding a scalar value if both vectors have the same number of elements. Appearing much more rarely, multiplying a column vector by a row vector is an outer product yielding a square matrix.

>

> Outer products have a common use case which isn't represented in the

> proposal: Rank-1 update of an existing matrix. If you try to spell

> that using overloaded arithmetic operators, you get the following,

> assuming that x is a row vector and y is a column vector.

(Did you mean x is a column and y a row?)

>

> A += x * y;

>

> A naive implementation would always create a new temporary matrix to

> hold the outer product result x * y. The only way NOT to do that, and

> still retain the syntax "A += x * y," would be to use expression

> templates. ("x * y" would return outer_product_expression<X, Y>, and

> matrix::operator+=(outer_product_expression<X, Y>&&) would perform a

> rank-1 update.)

Right. But expression templates arenâ€™t a panacea for matrix-based operations. E.g., if y is replaced by an expression that depends on A, you run into aliasing issues.

>

> For many users of my libraries, A is M x N where M and N are often

> much bigger than 4. Creating a temporary matrix on every rank-1

> update would be prohibitively expensive. BLAS 2 - style in-place LU

> factorization of an N x N matrix involves N - 1 outer products, and

> would therefore involve N - 1 allocations of max size (N-1)^2.

> However, LU factorization only needs kN extra space for a small

> positive integer k. Even for 4 x 4 matrices and length-4 vectors, it

> could save (more precious SIMD) registers and instructions not to

> create a temporary matrix to hold the outer product result.

>

> Therefore, if you want to retain use of operator* for outer products,

> I conclude that the library would need to promise to use expression

> templates. Otherwise, I conclude that the rank-1 update use case

> would be better served by a named function (e.g., outer(x, y, A)) that

> updates an existing matrix in place.

Agreed: I think a functional interface is more reliable at the outset.

Daveed

>

> Greetings and thanks for your hard work!

>

> I've made it a point to limit my comments about P1385. It's a

> different design that serves a different set of users. That being

> said, I'd like to comment on outer products and rank-1 updates.

>

> On Mon, May 1, 2023 at 5:55 AM Guy Davidson via SG19

> <sg19_at_[hidden] <mailto:sg19_at_[hidden]>> wrote:

>>

>> Hello everyone

>>

>> I'm just putting together a first pass at the wording for P1385, A proposal to add linear algebra support to the C++ standard library. If you look at the latest revision you will infer that at Kona in November and at Issaquah in February I addressed SG6 and LEWG about withdrawing the vector class entirely and simply offering a matrix class, where a vector is a special case of a matrix, with a single row or column. There were no objections to this approach.

>>

>> While there were no objections raised in the meeting, others have come in, and I want to use the reflectors to gather opinion about the matter. The heart of the problem is: what does the vector product signify? Is it an inner or outer product? Is vector orientation significant?

>>

>> With my mathematician's hat on, multiplying a row vector by a column vector is an inner product, yielding a scalar value if both vectors have the same number of elements. Appearing much more rarely, multiplying a column vector by a row vector is an outer product yielding a square matrix.

>

> Outer products have a common use case which isn't represented in the

> proposal: Rank-1 update of an existing matrix. If you try to spell

> that using overloaded arithmetic operators, you get the following,

> assuming that x is a row vector and y is a column vector.

(Did you mean x is a column and y a row?)

>

> A += x * y;

>

> A naive implementation would always create a new temporary matrix to

> hold the outer product result x * y. The only way NOT to do that, and

> still retain the syntax "A += x * y," would be to use expression

> templates. ("x * y" would return outer_product_expression<X, Y>, and

> matrix::operator+=(outer_product_expression<X, Y>&&) would perform a

> rank-1 update.)

Right. But expression templates arenâ€™t a panacea for matrix-based operations. E.g., if y is replaced by an expression that depends on A, you run into aliasing issues.

>

> For many users of my libraries, A is M x N where M and N are often

> much bigger than 4. Creating a temporary matrix on every rank-1

> update would be prohibitively expensive. BLAS 2 - style in-place LU

> factorization of an N x N matrix involves N - 1 outer products, and

> would therefore involve N - 1 allocations of max size (N-1)^2.

> However, LU factorization only needs kN extra space for a small

> positive integer k. Even for 4 x 4 matrices and length-4 vectors, it

> could save (more precious SIMD) registers and instructions not to

> create a temporary matrix to hold the outer product result.

>

> Therefore, if you want to retain use of operator* for outer products,

> I conclude that the library would need to promise to use expression

> templates. Otherwise, I conclude that the rank-1 update use case

> would be better served by a named function (e.g., outer(x, y, A)) that

> updates an existing matrix in place.

Agreed: I think a functional interface is more reliable at the outset.

Daveed

Received on 2023-05-01 21:43:29