[WIP] Add comparision of converged results to replication folder

[dmbates]

• Add script evaluation.jl to blockecholesky/replication/ with code for comparing evaluation of the objective for different processors, number of threads, and BLAS implementations.
• On Macs the PRIMA implementation of BOBYQA can be compared to the NLopt implementation. (It seems that this is not possible on amd_64 Linux implementations at present.)
• I think we should consider creating a section in the paper about how we are describing the evaluation of the objective but the actual fitting process involves numerical optimization of that objective in which minor differences can and do occur.

[dmbates]

@palday @ajinkya-k suggestions/modifications of the code and comparison format are welcome.

[dmbates]

I realized that it should be possible to start with OpenBLAS active then fit the model and extract information from it then load AppleAccelerate or MKL, as appropriate, then repeat the whole process.

@ajinkya-k If I recall correctly you had some code to check on the processor family and OS before loading an accelerated BLAS. Does this ring a bell?

[dmbates]

I added an evaluation summary from a linux_amd64 machine. It is not possible to use PRIMA on that computer. It seems that the PRIMA_jll package is not available.

[ajinkya-k]

@ajinkya-k If I recall correctly you had some code to check on the processor family and OS before loading an accelerated BLAS. Does this ring a bell?

Yes. I can dig it out

[dmbates]

@ajinkya-k If I recall correctly you had some code to check on the processor family and OS before loading an accelerated BLAS. Does this ring a bell?

Yes. I can dig it out

Thanks.

[ajinkya-k]

@dmbates Its there in the qmd already:

Manuscripts/blockedcholesky/BlockedCholeskyMM.qmd

Lines 454 to 460 in dba3bd5

	using MKL_jll
	hostarch = Base.BinaryPlatforms.arch(Base.BinaryPlatforms.HostPlatform())
	@static if Sys.isapple() && hostarch == "aarch64"
	using AppleAccelerate
	elseif MKL_jll.is_available()
	using MKL
	end

[dmbates]

@dmbates Its there in the qmd already:

Manuscripts/blockedcholesky/BlockedCholeskyMM.qmd

Lines 454 to 460 in dba3bd5

	using MKL_jll
	hostarch = Base.BinaryPlatforms.arch(Base.BinaryPlatforms.HostPlatform())
	@static if Sys.isapple() && hostarch == "aarch64"
	using AppleAccelerate
	elseif MKL_jll.is_available()
	using MKL
	end

Thanks.

[dmbates]

I added some discussion of replicability of the evaluation in a new section 3.4 I'm not sure to what extent we need to go over the fact that floating point arithmetic is not exact and one should expect differences from different BLAS implementations. I focussed on m.optsum.finitial because it isolates the difference in objective values to different BLAS implementations.

I hope this explanation doesn't result in recommendations to change the code to display fewer digits, which is the simple way to address problems of slightly different values on different machines. Well "simple" in concept, potentially quite difficult in implementation.

[ajinkya-k]

@dmbates sorry for the pedantic edit suggestion

[dmbates]

@ajinkya-k Not a problem. In programming pedantry is a good thing. 😄

[dmbates]

@ajinkya-k and @palday I added sections describing the unconstrained optimization and also an example of a likelihood ratio test. I haven't looked at the whole manuscript after these additions. There may be other sections that could be pared back now. Please take a look.

[ajinkya-k]

@dmbates I got this error:

ERROR: Julia server returned error after receiving "run" command:

Failed to run notebook: /Users/ajinkya/projects/julia_repos/Manuscripts/blockedcholesky/BlockedCholeskyMM.qmd

The underlying Julia error was:

EvaluationError: Encountered 1 error during evaluation

Error 1 of 1
@ /Users/ajinkya/projects/julia_repos/Manuscripts/blockedcholesky/BlockedCholeskyMM.qmd:612
ArgumentError: Log-likelihood must not be lower in models with more degrees of freedom
Stacktrace:
 [1] lrtest(::LinearMixedModel{Float64}, ::Vararg{LinearMixedModel{Float64}}; atol::Float64)
   @ StatsModels ~/.julia/packages/StatsModels/mPD8T/src/lrtest.jl:121
 [2] lrtest(::LinearMixedModel{Float64}, ::Vararg{LinearMixedModel{Float64}})
   @ StatsModels ~/.julia/packages/StatsModels/mPD8T/src/lrtest.jl:74
 [3] top-level scope
   @ ~/projects/julia_repos/Manuscripts/blockedcholesky/BlockedCholeskyMM.qmd:613


Stack trace:
    at writeJuliaCommand (file:///Applications/quarto/bin/quarto.js:78249:13)
    at eventLoopTick (ext:core/01_core.js:178:7)
    at async executeJulia (file:///Applications/quarto/bin/quarto.js:78126:20)
    at async Object.execute (file:///Applications/quarto/bin/quarto.js:77766:16)
    at async renderExecute (file:///Applications/quarto/bin/quarto.js:123336:25)
    at async renderFileInternal (file:///Applications/quarto/bin/quarto.js:123591:35)
    at async renderFiles (file:///Applications/quarto/bin/quarto.js:123385:9)
    at async render (file:///Applications/quarto/bin/quarto.js:129214:19)
    at async _Command.actionHandler (file:///Applications/quarto/bin/quarto.js:129456:24)
    at async _Command.execute (file:///Applications/quarto/bin/quarto.js:9833:7)

this comes from lrtest

[ajinkya-k]

i dont get this error when running the same code in the REPL

[dmbates]

I am using the current db/unconstrained version of MixedModels.jl for evaluating the code blocks, which may affect the convergence. Can you check if you are using that version? (I need to finalize the changes in MixedModels.jl on that branch and make a release.)

(blockedcholesky) pkg> st
Status `~/git/Manuscripts/blockedcholesky/Project.toml`
  [cbdf2221] AlgebraOfGraphics v0.11.3
  [13e28ba4] AppleAccelerate v0.4.1
  [69666777] Arrow v2.8.0
  [336ed68f] CSV v0.10.15
  [13f3f980] CairoMakie v0.15.5
  [0ca39b1e] Chairmarks v1.3.1
  [a93c6f00] DataFrames v1.7.0
  [33e6dc65] MKL v0.9.0
  [ff71e718] MixedModels v4.38.0 `~/.julia/dev/MixedModels`
  [0a7d04aa] PRIMA v0.2.4
  [08abe8d2] PrettyTables v2.4.0
  [856f044c] MKL_jll v2025.2.0+0
  [de0858da] Printf v1.11.0
  [2f01184e] SparseArrays v1.11.0

Also, I may have used the identifier m2 elsewhere so it could be important to clear any caches so that the m2 in the comparison is the m2 that was evaluated before that chunk. (I can't imagine how it could end up being anything else, but just in case ...)

[ajinkya-k]

I am not. but can that affect when rendering in qmd vs in repl?

[dmbates]

I don't know why there are different results in rendering the .qmd file versus execution in the REPL, other than the possibility that something was cached. I might try nuking the .cache directory.

Using an earlier version of MixedModels.jl, even the current released version, will use bounded optimization to fit the model, which may result in a different loglikelihood at convergence. I don't think it could be sufficiently different to cause this error but we might as well rule it out.

[ajinkya-k]

i tried cleaning the cache. removing all the generated files but the problem persists. I am using the release version of MixedModels though not the dev version

[dmbates]

@ajinkya-k I will email you the .pdf file that I generated so you can cross check the results from running the blocks in the REPL.

[dmbates]

Try turning off the evaluation of the block with the lrtest and seeing if the log-likelihoods correspond to those in the PDF file I sent. (The reason I suggest turning off evaluation of that block is just to get through to the end of the evaluation and production of the PDF.)

[ajinkya-k]

Will do!

[ajinkya-k]

the likelihood is different! (left is what i get, right is your version)

[dmbates]

In your version theta_3 goes to zero. Is it possible that you are using a version of MixedModels.jl with the non-negativity constraint? Here's the end of the fitlog for that model on my system

julia> m2.optsum.fitlog
490-element Vector{Tuple{Vector{Float64}, Float64}}:
 ([1.0, 1.0, 1.0, 0.0, 1.0], 242059.14387186497)
 ([2.0, 1.0, 1.0, 0.0, 1.0], 245732.8733592998)
 ([1.0, 2.0, 1.0, 0.0, 1.0], 243381.25249141201)
 ([1.0, 1.0, 2.0, 0.0, 1.0], 242078.25849461395)
 ([1.0, 1.0, 1.0, 1.0, 1.0], 242059.77726579554)
 ([1.0, 1.0, 1.0, 0.0, 2.0], 242078.26388572814)
 ([0.0, 1.0, 1.0, 0.0, 1.0], 241071.80487130157)
 ([1.0, 0.0, 1.0, 0.0, 1.0], 251890.34637374704)
 ([1.0, 1.0, 0.0, 0.0, 1.0], 242011.6120821771)
 ([1.0, 1.0, 1.0, -1.0, 1.0], 242059.26244893295)
 ([1.0, 1.0, 1.0, 0.0, 0.0], 242074.39366040094)
 ([0.1351951978318343, 1.3814470900881257, 0.0870282613136546, -0.006959733251804052, 0.9484777483061189], 239938.40680260953)
 ([-0.19110577769441361, 1.4816788474648885, -1.8762415514965225, -0.03145036060361427, 0.7797972234551669], 239579.96381637195)
 ([-0.2271537808848385, 1.4399358386530774, -2.977541151909495, -0.5626112738464235, -0.8019629786787126], 239349.85757374717)
 ([-0.21519531002297176, 1.4490839906869457, -4.801008006051292, -2.392674878368525, -3.855749012692688], 239468.17955458688)
 ([-0.3229264475029543, 1.4318141656070045, -4.185241374925724, 0.5744963958853866, -1.9151601550672597], 239330.6080274484)
 ([-0.2911669805392789, 1.4414591445521123, -3.7686986497139774, 1.4827531312990403, -1.9366117584626958], 239297.916398924)
 ⋮
 ([-0.27570356743452346, 0.4352620055644907, -0.046105475339223975, 0.07782694591791324, -0.11294087289572423], 237647.0585686748)
 ([-0.27573005035130355, 0.4352901076126246, -0.046113147967613725, 0.07783420678451904, -0.11284923589119428], 237647.0585265882)
 ([-0.2757031625959591, 0.43527418887654495, -0.046034864442280506, 0.07783159120749195, -0.11279549120892356], 237647.05845407228)
 ([-0.27569850182476885, 0.4352840731783276, -0.04601479325090905, 0.07779459396187309, -0.11270544161767668], 237647.05840824364)
 ([-0.2756916037258249, 0.43528245356901996, -0.04599645389069544, 0.07760375798035114, -0.11264891649830344], 237647.0584083805)
 ([-0.275639921399884, 0.43536461702920276, -0.04601094473105374, 0.07780266006610882, -0.11270436981363542], 237647.0585529526)
 ([-0.27569639166822313, 0.4352742329739572, -0.04599324063765421, 0.07777927796664454, -0.11260952695368129], 237647.05839018204)
 ([-0.27570283062700207, 0.43528028831534243, -0.04598301558258866, 0.07794544659567493, -0.11250870254275527], 237647.0583829815)
 ([-0.27569902828361786, 0.4352797224200382, -0.04590893515238784, 0.07800269578865285, -0.11247377631032274], 237647.05840374975)
 ([-0.2757111168189777, 0.43530858967567115, -0.045989820105484365, 0.07800172068299516, -0.11243177856996056], 237647.05839114104)
 ([-0.27567553568962766, 0.4353258387765646, -0.045937222956316685, 0.07797122516143218, -0.11244222967007185], 237647.05841263288)
 ([-0.27569570711449615, 0.4352884374851639, -0.04598857248599069, 0.07785152212937413, -0.11248133753748815], 237647.05838447565)
 ([-0.2757074813267906, 0.4353135297572702, -0.04593996427042295, 0.07802354585723645, -0.11253904174832151], 237647.05840307768)
 ([-0.2756973448847976, 0.43529192783313586, -0.04599680696609354, 0.07793861629066005, -0.11250127110474602], 237647.05837997136)
 ([-0.2756968028118065, 0.43529252198468454, -0.04600248894918995, 0.0779354916119583, -0.11250388233356183], 237647.05838002358)
 ([-0.27570596085578436, 0.4352962913153145, -0.04599564097179936, 0.07793633288574711, -0.11250166044590829], 237647.05838153942)
 ([-0.27569764077585507, 0.4352944305191494, -0.04599700916252852, 0.07793015589164233, -0.11250176488366956], 237647.0583798222)

[ajinkya-k]

That might be it. I might be using the latest stable version (v4.38) instead of the dev version

[dmbates]

I should have made a new release with the unconstrained optimization but I wanted to get a bit more experience with it first.

[dmbates]

Notice that model m2 is a case of converging outside the boundary area in that \theta_3 is negative. The signs of \theta_3 and \theta_4 will be flipped after convergence.

[ajinkya-k]

Suggested change

	MixedModels = "4.34"
	MixedModels = "4.38"

[ajinkya-k]

things worked after I updated the compat bounds and used the dev version. Of course we will have to update to 4.39 after that release is available