Matrix factorizations with unitful elements fail

[jiahao]

A simple example:

julia> using SIUnits; A=randn(1,1)Second
1x1 Array{SIQuantity{Float64,m,kg,s,A,K,mol,cd},2}:
 -1.0710545530502584 s 

julia> lufact(A)
ERROR: `convert` has no method matching convert(::Type{Float64}, ::SIQuantity{Float64,0,0,1,0,0,0,0})
 in copy! at abstractarray.jl:149
 in convert at abstractarray.jl:303
 in lufact at linalg/lu.jl:65

The error arises here because lufact dispatches on this method which computes T=eltype(A); S=one(T)/one(T) and converts the matrix to an AbstractMatrix{S}. In this case S == Float64 drops the SIQuantity annotation, which causes the error.

cc: @Keno @andreasnoack

[Jutho]

What would you expect for the factorisation of a unitful matrix. Where should the units go? In the L or in the U of the lufact? Or should they both carry, in this case, a unit s^(1/2)?

[jiahao]

It's not entirely obvious to me, but I'd propose for a matrix M with units s:

• for ordinary LL^T Cholesky, L should have units s^(1/2)
• for LDL^T Cholesky, L should be unitless and D has units s
• for LU, both L and U have units s^(1/2)
• for QR, Q is unitless and R has units s

[andreasnoack]

The generic factorization function almost works for a matrix with seconds. As a natural consequence of the algorithm U will be in seconds and L will not have units. There are two one problems:

• There is no real(SIQuantity) method, but that should be easy enough to define.
• (Wrong. Changed my mind on this one) inv(x::Number) doesn't work correctly for SUQuantitys, instead of one(x)/x I think it should be 1/x.

When I fix the two problem s, the factorization works. Then there is the promotion problem. For the factorization to work when the elements have units, it is necessary that the array is heterogenous, which is also the case when you define it by

A=randn(3,3)Second;

but the problem is that one(T)/one(T) drops the type. The promotion problem could either be solved by using the ugly rule promote_type(typeof(one(eltype(A))/one(eltype(A))), eltype(A)) or by not dropping the SIQuantity type when dividing.

[andreasnoack]

It is a bit tricky with inv(x::Number). Actually, I think the definition is correct because one(x) should be the multiplicative identity. However, this is not the case for one(x::SIQuantity) because

julia> x=1.0Second
1.0 s 

julia> x*one(x)
1.0 s²

julia> x==x*one(x)
false

[StefanKarpinski]

That's a tough one.

[Jutho]

Even for an x::SIQuantity, one(x) should have no unit. It should indeed be the multiplicative identity. Something with units cannot be an identity, as @andreasnoack clearly illustrated above. Another way to see this is that an identity should not depend on the chosen units, and clearly 1s = 1000 ms cannot be an identity.

[StefanKarpinski]

Of course zero(si) needs to still have units.

[alanedelman]

In LU , L should be unitless and U should have units s
Follow the pivots
I think

On Tue, Jul 15, 2014 at 11:06 PM, Jiahao Chen [email protected]
wrote:

It's not entirely obvious to me, but I'd propose for a matrix M with units
s:

• for ordinary LL^T Cholesky, L should have units s^(1/2)
• for LDL^T Cholesky, L should be unitless and D has units s
• for LU, both L and U have units s^(1/2)
• for QR, Q is unitless and R has units s

—
Reply to this email directly or view it on GitHub
#128.

[stevengj]

@jiahao, note that SIQuantity currently only allows each SI unit to have an integer power, so fractional powers aren't currently possible. This makes units for Cholesky factorizations problematic.

[stevengj]

I agree with Alan about A=LU. (As @andreasnoack pointed out, Gaussian elimination produces a dimensionless L.)

[andreasnoack]

I believe this is now just an SIUnits issue. It has

julia> A = eye(2)Second                                           
2x2 Array{SIUnits.SIQuantity{Float64,m,kg,s,A,K,mol,cd,rad,sr},2}:
 1.0 s  0.0 s                                                     
 0.0 s  1.0 s                                                     

julia> typeof(one(A[1]))
Float64                                                           

but

julia> typeof(one(eltype(A)))                
SIUnits.SIQuantity{Float64,0,0,0,0,0,0,0,0,0}

because it falls back on the Number fallback. It's a instance of a type for which the right thing is typeof(x) != typeof(one(x)). Char is similar but for zero since

julia> 'a' == 'a' + 0
true

(but we don't define zero(Char)).

Right now, qrfact doesn't work, but that is because of a missing abs2 in SIUinits.

cc: @Keno