Proposed LWG issue for `vector_sum_of_squares` and related functions

[MattStephanson]

I think there are a few issues with the vector_sum_of_squares wording, in addition to #141 and #144. The essence is that (1) the exact value of the scaling factor ought to be QoI detail and (2) there's no formal requirement that division or square root be defined for the Scalar type.

I've drafted the following LWG issue. I was wondering if anyone has any comments before I submit it, or would like to collaborate on writing it.

---

Discussion

The Effects of vector_sum_of_squares are both too strong and too weak. They are too weak because that function is not a triangular solve algorithm, so there is no requirement that it be possible to divide two values, and without division it is impossible to calculate the scaled sum-of-squares. They are too strong because they describe precisely what value result.scaling_factor must have, but implementations should have freedom to choose this as needed to avoid unnecessary under- or overflow. Notably, the reference BLAS stopped using this prescription in 2021 (see Reference-LAPACK/lapack#494). An implementation may also prefer to scale by a power of the floating-point radix, because this can be done exactly.

The specifications of [linalg.algs.blas1.nrm2] and [linalg.algs.blas1.matfrobnorm] are also too weak because they could mandate that it be possible to take the square root of a Scalar but don't, and that operation is inherent in the function's definition.

Proposed resolution

[linalg.general]

-7- Within all the functions in [linalg], square roots shall be calculated with an unqualified call to sqrt, with overload resolution performed in a context that includes the declaration template<class U> U sqrt(const U&) = delete;.

[linalg.reqs.alg]/2.2

The triangular solve algorithms in [linalg.algs.blas1.ssq], [linalg.algs.blas2.trsv], [linalg.algs.blas3.trmm], [linalg.algs.blas3.trsm], and [linalg.algs.blas3.inplacetrsm] either [...]

[linalg.algs.blas1.ssq]/3

Effects: Returns a value result such that

(3.1) — result.scaling_factor is the maximum of init.scaling_factor and abs-if-needed(x[i]) for all i in the domain of v; and

(3.2) — letif s2init be is

init.scaling_factor * init.scaling_factor * init.scaled_sum_of_squares

then [...]

[linalg.algs.blas1.nrm2]/2

Mandates: decltype(sqrt(abs-if-needed(declval<typename InVec​::​value_type>()))) is convertible to Scalar.

[linalg.algs.blas1.matfrobnorm]

Mandates: Let a be abs-if-needed(declval<typename InMat​::​value_type>()). Then, decltype(sqrt(init + a * a)) is convertible to Scalar.

[mhoemmen]

Thanks for pointing this out! I think we meant to file an LWG issue, but it got dropped among all the other things to do.

The Effects of vector_sum_of_squares are both too strong and too weak. They are too weak because that function is not a triangular solve algorithm, so there is no requirement that it be possible to divide two values, and without division it is impossible to calculate the scaled sum-of-squares.

The original intent was for vector_sum_of_squares to behave like LAPACK's xLASSQ for floating-point types, and otherwise to compute the sum of squares naively and set scaling_factor to one. However, the current wording breaks this in two ways.

1. Requiring a particular value of scaling_factor requires more operations, like sqrt and division, as you point out.
2. "Set[ting] scaling_factor to one" requires a way to get the multiplicative identity for Scalar.

Regarding Breakage (1), we want vector_sum_of_squares to work for all types, including types that don't have sqrt or the division operator. For example, for an mdspan of int, it shouldn't try to take square roots or do division. The problem is that argument-dependent lookup (ADL) would find sqrt for an int argument. As a result, an implementation of your suggested wording would break vector_sum_of_squares for mdspan of int. Likewise, division compiles for int but would do the wrong thing in vector_sum_of_squares.

Regarding Breakage (2), the only default value that would satisfy [linalg.algs.blas1.ssq] 3.2 in general is the multiplicative identity. Imposing the triangular solve algorithms requirements (per [linalg.algs.req] 2.2) would not help here. In general, std::linalg offers no way to construct the multiplicative identity. Fixing this would require a design change.

Given that the C++26 design deadline has passed, the least intrusive way to fix this issue would be to Mandate InVec​::​value_type and Scalar to be arithmetic types or specializations of complex. It's never a breaking change to go from "ill-formed if X" to "well-formed and well-defined behavior if X." We would then propose a design fix for C++29 that would make vector_sum_of_squares have the desired behavior for more general types.

Would that be an acceptable change?

[mhoemmen]

Notably, the reference BLAS stopped using this prescription in 2021 (see Reference-LAPACK/lapack#494). An implementation may also prefer to scale by a power of the floating-point radix, because this can be done exactly.

I missed this point in my discussion above. This is a good reason to pull the existing wording or even pull the entire function.

[MattStephanson]

[Sorry, I drafted this yesterday, before your second message, and forgot to post it.]

Hi Mark,

1. Limiting vector_sum_of_squares to arithmetic types or specializations of complex sounds okay to me, though if it could be done in a way that also permits user-defined complex types that would be better. It would be convenient of there were already an exposition-only concept for that.
2. [Per your second comment above,] I do think it's important to address the overspecification of scaling_factor sooner, if we want to address it at all. If it gets delayed to C++29, then it might be much harder to relax later. I understand the reference to xLASSQ and its earlier behavior, but as I pointed out using the maximum value of the inputs hasn't been used even in the reference BLAS for several years.
3. My sqrt comments were specifically for vector_two_norm and matrix_frob_norm, not vector_sum_of_squares. The former two specifically mention "square root" in their Returns clauses.

[mhoemmen]

Greetings @MattStephanson ! Just wanted to give you an update.

I proposed the vector_sum_of_squares issue to LWG and it was filed as LWG 4302. Given that it's a design change, it's likely that the issue would be referred back to LEWG. I plan to submit an NB comment on this. If the NB comment is not accepted, the best I can do is advocate that the design change be accepted retroactively for C++26 as a Defect Report (DR).

I will file another LWG issue for vector_two_norm and matrix_frob_norm. As you pointed out, the wording is underspecified, since it does not explain how to take a "square root." I'll post the LWG issue number here once it's up.

[MattStephanson]

Thanks @mhoemmen !

[mhoemmen]

Just FYI, here is the proposed issue I filed. Thanks for your input! : - )

Issues

The wording for vector_two_norm [linalg.algs.blas1.nrm2] and matrix_frob_norm [linalg.algs.blas1.matfrobnorm] has two issues.

1. Their Returns clauses say that the functions return the "square root" of the sum of squares of the initial value and the absolute values of the elements of the input mdspan. However, nowhere in [linalg] explains how to compute a square root.

   1.a. The input mdspan's value_type and the initial value type are not constrained in a way that would ensure that calling std::sqrt on this expression would be well-formed.

   1.b. There is no provision to find sqrt via argument-dependent lookup, even though [linalg] has provisions to find abs, conj, real, and imag via argument-dependent lookup. There is no "sqrt-if-needed" analog to abs-if-needed, conj-if-needed, real-if-needed, and imag-if-needed.

2. The overloads that take an initial value parameter Scalar init return Scalar.

   2.a. This may silently lose information if the function uses std::sqrt to compute square roots. For example, if Scalar and the input mdspan's value_type are both int, the square root computed via std::sqrt would return double. However, vector_two_norm and matrix_frob_norm returning Scalar would force a rounding conversion back to int.

Suggested fix

The easiest fix for both issues is just to Constrain both Scalar and the input mdspan's value_type to be floating-point numbers or specializations of std::complex for these two functions. This presumes that relaxing this Constraint and fixing the above two issues later would be a non-breaking change. If that is not the case, then I would suggest removing the two functions entirely.

Related non-issues

The return type being the same as the initial value parameter type (Scalar) has the following additional implications. We do not consider these wording bugs. We state them here to show future readers that we did consider these implications.

In what follows, let a "user-defined complex number type" be a linear algebra value type T that is neither an arithmetic type nor a std::complex specialization, such that

• abs with argument declval<T>() is findable via argument-dependent lookup, and

• the type of abs(declval<T>()) differs from remove_cvref_t<T>.

2.b. If Scalar is a specialization of std::complex or a user-defined complex number type, then vector_two_norm and matrix_frob_norm would return a complex number. Even though that number would have zero imaginary part, these functions mathematically should return magnitudes, so their return type should be non-complex.

2.c. If Scalar is a user-defined complex number type, then (2.b) implies that the return type of abs must be convertible to Scalar. However, [linalg] does not constrain user-defined complex number types to make this conversion well-formed.

We do not consider these issues because the initial value mathematically represents a previous partial result -- a magnitude. It should not be a complex number type. Users are free to use std::complex<Real> to represent magnitudes; the result would be well-formed, its real part would be the desired magnitude, and its imaginary part would be zero (due to conversion from Real to std::complex<Real>). If a user-defined complex number type doesn't have a conversion from a magnitude, then users shouldn't be using it as the initial value type. The same argument applies to other functions in [linalg] that return absolute values or sums of absolute values, such as vector_sum_of_squares, vector_abs_sum, matrix_one_norm, and matrix_inf_norm. (For vector_sum_of_squares, please note LWG Issue 4302.)

Neither (2.b) nor (2.c) apply to the overloads that do not take an initial value parameter, because they deduce the return type as the result of abs-if-needed(declval<value_type>()).

Suggested wording change

1. In [linalg.algs.blas1.nrm2] (Euclidean norm of a vector), insert the following text as a new paragraph between 1 and 2: Constraints: InVec::value_type and Scalar are either a floating-point type, or a specialization of std::complex.

2. In [linalg.algs.blas1.matfrobnorm] (Frobenius norm of a matrix), insert the following text as a new paragraph between 1 and 2: Constraints: InVec::value_type and Scalar are either a floating-point type, or a specialization of std::complex.

Thanks to Matt Stephanson for pointing out Issue (1) (see #308 for discussion).