Too restrictive setindex! for triangular matrices

[mateuszbaran]

setindex! of triangular matrices could be less restrictive when setting to a structural zero: https://github.com/JuliaLang/julia/blob/4ed484fe7ac591da67bfe384c369430ed658606d/stdlib/LinearAlgebra/src/triangular.jl#L231
it could be iszero(x) instead.

I'm not sure about diagonals of unit triangular matrices (Julia doesn't have isoneunit) but isone(x) may be better then x == 1 as well.

[stevengj]

iszero(x) here seems fine, but for any numeric type x == 0 and iszero(x) should be functionally equivalent, no?

[mateuszbaran]

Yes, I haven't encountered any numeric types for which these two are different. What I was doing is making tests for JuliaArrays/StaticArrays.jl#814 using block matrices and I was surprised that LinearAlgebra didn't work for block-triangular matrices. If you think triangular matrices are only useful for numeric types, feel free to close this issue 🙂 .

[stevengj]

Triangular matrices are certainly useful in principle for block-triangular types, but I suspect that a number of code updates are required to make that work with our existing types. Using iszero here is a good start, of course, and is certainly harmless at worst.

[mbauman]

We've had this same discussion in the context of SparseArrays; we should match that. JuliaLang/julia#30580

[mateuszbaran]

We've had this same discussion in the context of SparseArrays; we should match that. JuliaLang/julia#30580

All that PR does is changing == 0 to iszero where appropriate. This could be easily extended to UpperTriangular and LowerTriangular (and UpperHessenberg as well). I'm not sure what should be done about diagonal of unit-triangular matrices to make it consistent. Perhaps == oneunit(T)? cc @andreasnoack (we had a short discussion on Slack about one vs oneunit for unit-triangular matrices).

[stevengj]

x == oneunit(x) is the correct comparison, and will work if x is a matrix. (However, it will be suboptimal compared to a hypothetical isone function since oneunit(x) will instantiate a new matrix for comparison.) oneunit(T) won't work for matrix types, because it doesn't know the size of the matrix. oneunit, not one, is more correct here since matrices can be dimensionful.

So I would support changing x == 0 to iszero(x) and x == 1 to x == oneunit(x). It shouldn't hurt anything in the numerical-entry case.

[andreasnoack]

Using iszero here should be uncontroversial. However, I disagree that oneunit is the right choice here. Possibly the most common unit triangular matrix is the L matrix in LU and we used to agree that U should carry the units and L should be unitless.

[mateuszbaran]

Currently getindex for unit-triangular matrices returns oneunit(T) for the diagonal: https://github.com/JuliaLang/julia/blob/2a961764962f6044416c1e99d012785066333ac9/stdlib/LinearAlgebra/src/triangular.jl#L220-L221
It also returns zero(T) outside the diagonal which will fail on many matrix types (for block-triangular matrices) as @stevengj pointed out. Just guessing but maybe something like zero(transpose(A.data[j,i])) would be fine? Maybe LinearAlgebra should offer both variants of unit-triangular matrices, that is unitful and unitless diagonal?

[stevengj]

Using iszero here should be uncontroversial. However, I disagree that oneunit is the right choice here. Possibly the most common unit triangular matrix is the L matrix in LU and we used to agree that U should carry the units and L should be unitless.

If x is dimensionless then oneunit(x) is equivalent to one(x).

The point is, whatever units L has (if any), you should be calling setindex! with a value in the same units (if any), in which case oneunit is the right thing. Though maybe it should be oneunit(L.data[i,i])

(Note that iszero(x) also uses the same units as x to compare with zero, if x is dimensionful.)