rankUpdate for BlockDiagonal by Dense #447
Conversation
Codecov Report
@@ Coverage Diff @@
## master #447 +/- ##
==========================================
+ Coverage 92.74% 92.85% +0.11%
==========================================
Files 23 23
Lines 1722 1736 +14
==========================================
+ Hits 1597 1612 +15
+ Misses 125 124 -1
Continue to review full report at Codecov.
|
|
@kleinschmidt While working on this, I discovered a second issue. These two formulae should give the same model: julia> fnested = @formula(Thickness ~ 1 + Source + (1+Source|Lot/Wafer))
FormulaTerm
Response:
Thickness(unknown)
Predictors:
1
Source(unknown)
(Source,Lot,Wafer)->(1 + Source) | Lot / Wafer
julia> f = @formula(Thickness ~ 1 + Source + (1+Source|Lot) + (1+Source|Lot&Wafer))
FormulaTerm
Response:
Thickness(unknown)
Predictors:
1
Source(unknown)
(Source,Lot)->(1 + Source) | Lot
(Source,Lot,Wafer)->(1 + Source) | Lot & WaferBut julia> fit(MixedModel, fnested, MixedModels.dataset(:oxide))
ERROR: MethodError: no method matching +(::Tuple{InterceptTerm{true},CategoricalTerm{DummyCoding,String,1}}, ::Tuple{RandomEffectsTerm,RandomEffectsTerm})
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
+(::MutableArithmetics.Zero, ::Any) at /home/XXX/.julia/packages/MutableArithmetics/0tlz5/src/rewrite.jl:52
+(::Tuple{Vararg{Union{AbstractTerm, Tuple{Vararg{AbstractTerm,N}} where N},N}} where N, ::AbstractTerm) at /home/XXX/.julia/packages/StatsModels/cpRh6/src/terms.jl:404
...
Stacktrace:
[1] +(::InterceptTerm{true}, ::CategoricalTerm{DummyCoding,String,1}, ::Tuple{RandomEffectsTerm,RandomEffectsTerm}) at ./operators.jl:538
[2] sum(::Tuple{InterceptTerm{true},CategoricalTerm{DummyCoding,String,1},Tuple{RandomEffectsTerm,RandomEffectsTerm}}) at ./tuple.jl:396
[3] apply_schema(::Tuple{ConstantTerm{Int64},Term,FunctionTerm{typeof(|),var"#107#109",(:Source, :Lot, :Wafer)}}, ::MixedModels.MultiSchema{StatsModels.FullRank}, ::Type{T} where T) at /home/XXX/Work/MixedModels.jl/src/schema.jl:30
[4] apply_schema(::FormulaTerm{Term,Tuple{ConstantTerm{Int64},Term,FunctionTerm{typeof(|),var"#107#109",(:Source, :Lot, :Wafer)}}}, ::StatsModels.Schema, ::Type{LinearMixedModel}) at /home/XXX/Work/MixedModels.jl/src/schema.jl:53
[5] LinearMixedModel(::FormulaTerm{Term,Tuple{ConstantTerm{Int64},Term,FunctionTerm{typeof(|),var"#107#109",(:Source, :Lot, :Wafer)}}}, ::NamedTuple{(:Source, :Lot, :Wafer, :Site, :Thickness),Tuple{Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.Primitive{Float64,Array{Float64,1}}}}; contrasts::Dict{Symbol,Any}, wts::Array{Any,1}) at /home/XXX/Work/MixedModels.jl/src/linearmixedmodel.jl:67 |
src/linalg/rankUpdate.jl
Outdated
| for i in 1:blksize, j in 1:blksize | ||
| iind = ioffset + i | ||
| jind = joffset + j | ||
| Arow = view(A, iind, :) |
There was a problem hiding this comment.
I would probably write this as an explicit loop over axes(A, 2) and avoid the construction of views.
There was a problem hiding this comment.
So something like this:
AtAij = 0
for idx in axes(A, 2)
# because the second is actually A', we swap index orders
AtAij += A[iind,idx] * A[jind,idx]
end
Cdat[i,j,k] += α * AtAij?
I also thought about writing this as
AtAij = sum(A[iind,idx] * A[jind,idx] for idx in axes(A,2))but then I realized the construction of the generator is probably no better (and maybe worse) than the construction of views.
There was a problem hiding this comment.
Actually what I had in mind is something like
for idx in axes(A, 2)
Cdat[i,j,k] += α * A[iind, idx] * A[jind, idx]
end
I'm pretty sure the compiler will detect that the location Cdat[i,j,k] is a loop invariant.
There was a problem hiding this comment.
Ah, interesting.
That might be a decent research project if we ever get an intern -- taking a look at @code_llvm and @code_native for some of our loops and seeing where the Julia compiler and LLVM are doing nice vectorization and other loop transformations for.
There was a problem hiding this comment.
An alternative is to resize Cdat to be a matrix with blksize rows. Then the columns match up with the columns of C. Also there is an intriguing function called fldmod1 which makes it easy to step over the blocks but I have forgotten the details. I'll take a look.
There was a problem hiding this comment.
If we're going to resize in that way -- would it make sense to change the internal storage used by UniformBlockDiagonal to always store in that format?
There was a problem hiding this comment.
There are times that you want the diagonal blocks to be accessed as k by k matrices. The most important part is in evaluating Lambda'Z'ZLambda + I. Also, the solving of multiple triangular systems to produce the [i,1] blocks of L for i > 1. I think it is easier to resize from 3 dimensions to 2 when we need it, rather than the other way around.
Co-authored-by: Douglas Bates <[email protected]>
Co-authored-by: Douglas Bates <[email protected]>
Co-authored-by: Douglas Bates <[email protected]>
Co-authored-by: Douglas Bates <[email protected]>
|
That's really strange, I don't see why that wouldn't work. It looks like it's not failing at the the random slope point, but rather at the point where it's trying to combine the two RE terms with the FE terms...but that path should already be hit by the existing tests, so I'm confused. I'll try to take a look when I get a chance. |
|
@kleinschmidt Thanks! Glad to know at least the problem isn't something obvious I should have been able to fix in 10 minutes... |
Another one from #445.