Support committee models

[davkovacs]

It would be nice to have a coefficient matrix instead of vector, where energy uses the first column (the mean that we use now), but there are further N columns that store committee parametrisations.

This needs some changes in IPFitting :brr solver as well.

[cortner]

what about a vector of svectors

[casv2]

Thanks David, I'm unsure how much needs to change. We now save/read an IPSuperBasis, but perhaps it'd be be able to store and evaluate the basis only and dot in the mean/committees. This will require more changes, and I'm not sure how deep we want to go. Maybe there's a clever hack?

[cortner]

and how will you then implement the evaluate function? will it return a vector, or a number with uncertainty, or ???

[casv2]

Exactly it's not obvious

[casv2]

We may have to experiment with this, we can take the mean norm differences of the committee members relative to the mean, that seems to work well in HAL

[cortner]

mutable struct PIPotential{T, NZ, TPI, TEV} <: SitePotential
   pibasis::TPI
   coeffs::NTuple{NZ, Vector{T}}
   dags::NTuple{NZ, CorrEvalGraph{T, Int}}
   evaluator::TEV
end

See - it's really easy to just have T = SVector{NC, Float64}. But what you do with it is another question. So I think we should flesh this out before diving in.

[davkovacs]

I would say it should evaluate to a number (the mean) or a vector of numbers if the committee members are queried.

[cortner]

so maybe evaluate should just give you the mean, and if you want something else you it should be done with a different function? I don't know; just a thought. But OTOH, if you want just the mean, then you don't need this. So maybe the committee potential should always return the error as well.

[casv2]

Yeah in principle returning the vector of numbers is fine for now at least, it'll actually allow people to explore it a bit more and maybe run some dynamics and monitor the committee members and see how they behave.

[davkovacs]

Exactly, but you wouldn't use anything else but the real mean for dynamics.

[casv2]

Someone can make a few nice plots, maybe discuss the range of different things we can do with them, discuss them in a group meeting and set on something that we all like

[davkovacs]

Next software meeting agenda?

[cortner]

we can have a WIP-PR ready quickly. One concern I have is how the action of evaluate propagates through the JuLIP energy and forces - this is all non-trivial I'm afraid.

[LarsSchaaf]

In addition to the N specific committees, wouldn't it be be nice if we just save the covariance matrix of the weights that is supplied by some solvers (like BRR) out of the box. So we can determine how many committees we want later on. I guess this would make the .json file quite large - which is why this shouldn't be the default behaviour?

[casv2]

Committees were never the plan, that's why, I'm also not convinced they have a particularly useful applications apart from usage in some database generation scheme to be honest. Adds cost and gives a reasonably good indication of error but it's far from perfect.

Why would we want to store the covariance matrix again? To evaluate the uncertainty in the full posterior predictive or to sample from? It's shaped Nbasis x Nbasis so pretty massive and uncomfortable to store. What's the benefit over just having a few committee members for people to try out?

[casv2]

so maybe evaluate should just give you the mean, and if you want something else you it should be done with a different function? I don't know; just a thought. But OTOH, if you want just the mean, then you don't need this. So maybe the committee potential should always return the error as well.

I think we should make efforts to do this, with the struct mentioned above and a function just returning the committee members for now, people can then do whatever to try things out?

[cortner]

does anybody want to try this: https://github.com/JuliaPhysics/Measurements.jl

[cortner]

if this could be propagated through everything, then we'd suddenly have error bars attached to everything we compute.

[cortner]

Just to record what Cas and I just discussed - it seems best to just add a committee parameter to the PIpot structure and implement a new evaluation code to make this efficient.

[cortner]

apologies this has been dormant. I'll think about this in the coming weeks though.

[cortner]

Ok - let's finally wrap this up. @casv2 and @MatthiasSachs and @WillBaldwin0 I need to know from you how you want the interface for this. Independently whether we use PIPot with a committee parameter or we define a new type LinearCommitteePotential, then I think evaluate, energy, forces, virial etc should behave exactly as before, but we can provide some additional functionality which does the same + uncertainties? E.g.

u_evaluate
u_evaluate_d
u_energy
u_forces
u_virial

And what should that return? the committee? The mean and variance? ... Do you want to call it something else?

Please specify exactly what you need and I'll then implement it. If we are sufficiently unclear on what we need then I think I will not put this in ACE1.jl after all. But let's see.

[casv2]

I think we need some ACE1.jl functions like com_energies(), com_forces(), com_virials() which just return K committee predictions in an array. Leveraging the 'dot-product trick' this should in principle incur (almost) no extra computational costs. It'd be nice to save/load ACE1.jl models with a T = SVector{K, Float64} containing the parameterisations too.

This requires minimal changes to ACE1.jl I think. We can compute means and variances to do UQ in some other package elsewhere. Same goes for setting up the biasing force which is an energy weighted force, requiring the committee predictions for both the energy and force. The (relative) force stopping criterion too, as well as other (new) biasing terms don't need to be computed inside ACE1.jl.

Just having access to the K predictions from ACE1.jl keeps things very modular and allows us to try out new ideas and leverage them to do UQ/HAL. I think this also means we wouldn't need any further changes inside ACE1.jl.

All of this is ofcourse would be besides the standard energy(), forces() stuff we have already corresponding to mean parameterisation/prediction provided by the BRR fit (maximum likelihood solution). It'd be nice to get both the mean (standard forces()) and K committee predictions (com_forces()) and only compute the basis evaluation once to make HAL really fast

So... maybe com_forces() should return K + 1 predictions, with the 1 index being the mean/maximum likelihood (or BRR fit) solution. I'd like to have access to the mean and committee predictions using a single and rapid call with almost no extra computational costs

[MatthiasSachs]

I agree with Cas. If we don't allow this flexibility, we would need specific functions for expectations of non-linear transformation, e.g., Expectation of ||Force on i-th atom||^2 and Expectation of ||Force on i-th atom|| or other statistics such as cov(complete force vector) . My only concern is memory efficiency in the case of large committees...

@casv2 are you currently using site-specific uncertainties? If so, I guess we may also want something like com_site_energies() which would return a vector of length #atoms with each entry being a committee of the respective site energy predictions?

[casv2]

Memory efficiency is a concern indeed, especially when having large systems and large committee members. However generally speaking I don't think we'll go much beyond say 32 members. I don't think we need to converge the Monte Carlo estimate much further as generally speaking the UQ isn't that accurate anyways...

Yes ofcourse com_site_energies() is much more relevant than com_energy() for UQ!

[cortner]

The problem with returning the committee is that you then can't exploit backpropagation. For energies this doesn't matter but for forces it matter very much.

[casv2]

Hmm... maybe we can leverage some threading? This is what I currently do in my version actually to speed things up

In practice what you'll see is that you end up running these HAL simulations on big nodes anyways with plenty empty threads. It works pretty well actually, the MD/HAL steps are pretty fast relative to say a DFT evaluation.

If you don't think this can be done faster by modifying ACE1.jl we may need to keep managing the committees ourselves in the different packages?

[cortner]

Sure but - we can always hack it - that’s not really the point.

If you take the N committees and then do a nonlinear transformation on them that spits out a scalar then you can backpropagate and get the gradient at the same cost as the evaluation. And this would probably even apply to forces - but if need to think this through first...

[cortner]

That said - a first implementation could be more limited. So let me ask this. To run HAL right now with the specific biases that work best what is the minimum you need?

[casv2]

Well so currently I use threading and manage my own committee essentially to set up the biasing. We can keep doing this for now because it works OK I guess? In that case we don't need any changes inside ACE1.jl for now...

We can just set up manage our committees in the agnostic HAL package (parameterised using ACEfit) and think of different biasing schemes in that package. Those can then be leveraged in some ACEhal package which performs the MD/MC steps?

It'd be nice to save the potentials with the committees, but it can wait perhaps until you've figured out the backprop with the force evaluation? It'd be very nice to have this though because HAL is then properly free and we can have all kinds of fun with very fast biased simulations

[cortner]

The agnostic HAL will only get the bias of force and or potential. It doesn't even need to know about committees I think?

I feel this needs another zoom call Monday or Tuesday?

[cortner]

cf #40

[cortner]

closed via #40 (backpropagation missing, making it a separate issue)