pathfinder missing lbfgs_status

[velochy]

In running pathfinder with a more complex generated model, I get the following error:

File ~/miniconda3/envs/salk/lib/python3.12/site-packages/pymc_extras/inference/pathfinder/lbfgs.py:190, in LBFGS.minimize(self, x0)
    [188](https://file+.vscode-resource.vscode-cdn.net/home/velochy/salk/sandbox/debug/kov1-full_fix_package/~/miniconda3/envs/salk/lib/python3.12/site-packages/pymc_extras/inference/pathfinder/lbfgs.py:188)     lbfgs_status = LBFGSStatus.CONVERGED
--> [190](https://file+.vscode-resource.vscode-cdn.net/home/velochy/salk/sandbox/debug/kov1-full_fix_package/~/miniconda3/envs/salk/lib/python3.12/site-packages/pymc_extras/inference/pathfinder/lbfgs.py:190) return history.x, history.g, history.count, lbfgs_status

UnboundLocalError: cannot access local variable 'lbfgs_status' where it is not associated with a value

Looking at the code in lbfgs.py

        # warnings and suggestions for LBFGSStatus are displayed at the end
        if result.status == 1:
            lbfgs_status = LBFGSStatus.MAX_ITER_REACHED
        elif (result.status == 2) or (history.count <= 1):
            if result.nit <= 1:
                lbfgs_status = LBFGSStatus.INIT_FAILED
            elif result.fun == np.inf:
                lbfgs_status = LBFGSStatus.DIVERGED
        else:
            lbfgs_status = LBFGSStatus.CONVERGED

        return history.x, history.g, history.count, lbfgs_status

There is a clear gap in the status==2 block that is causing this, but I don't understand the lbfgs algorithm nearly well enough to know what status that gap should be filled with.

[aphc14]

Hi @velochy, thanks for raising this issue. I'll try to address this one quickly. Can you set a breakpoint on this line:

elif (result.status == 2) or (history.count <= 1):

and share what output you get for result and history.count?

[velochy]

I'm sorry, but I have managed to misplace the model I actually got this error on :(

I do have cases of "All paths failed." and "BUG: Failed to iterate!" if you are interested.
Although they are generated models so it is a bit of work to get them working as stand-alone models, I'd be happy to try to get minimal examples if you (or someone else) are interested in taking them on
(before you ask - these are models that NUTS samples with relative ease, so the models themselves should be fine)

[aphc14]

Hi @velochy , would the minimal examples be too much hassle? If it's a bit of work to provide something, we could cover the gap with an else statement. With a minimal example, though, I can better understand:

1. If there should be 1 or more elif statements that can be added before an else, while letting the user know specifically what the issue is related to that individual Path, or if it's an issue related to the Paths.
2. If there needs to be another condition or filter of some sort when storing each point along the L-BFGS path
3. Types of scenarios where PF fails miserably, while NUTS can easily provide reasonable estimates (assuming it still results in zero successful paths after covering the gaps with an else statement).

All good if its too much work to create a minimal example that provides that error -- it's always possible to cover the gap now, and make some more changes later on.

[velochy]

I tried to extract it, and I think I solved my own problem.
I had copy-pasted code from an example that set jitter to 12... which clearly seems to be too large of a value for my data. Reducing it down to the default made the models sample nicely.

Anyways, maybe it is helpful for you - because I managed to reproduce both "BUG: Failed to iterate" as well as the lbfgs status errors with different random seeds (2 for lbfgs, 4 for failed to iterate)

import numpy as np, pymc as pm, pytensor.tensor as pt, pymc_extras as pme
inp = np.array([0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 2, 0, 0, 1, 1, 0, 1, 1,
       0, 1, 1, 0, 1, 0, 2, 0, 1, 0, 0, 0, 0, 1, 1, 1, 2, 0, 1, 2, 1, 0,
       1, 0, 1, 0, 1, 0])

res=np.array([[0,0,1,0,0],[0,0,1,0,0],[0,0,1,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,1,0,0,0],[0,0,0,1,0],[0,0,0,1,0],
[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[1,0,0,0,0],[0,0,1,0,0],[0,1,0,0,0],[0,0,0,1,0],[0,0,1,0,0],[0,0,1,0,0],
[0,0,0,1,0],[0,0,0,1,0],[0,1,0,0,0],[0,0,0,1,0],[0,0,1,0,0],[0,1,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[0,0,1,0,0],
[0,0,1,0,0],[0,0,1,0,0],[0,0,1,0,0],[0,1,0,0,0],[0,1,0,0,0],[0,0,0,1,0],[0,0,1,0,0],[0,1,0,0,0],[0,0,1,0,0],
[0,0,0,1,0],[0,0,1,0,0],[1,0,0,0,0],[1,0,0,0,0],[0,0,1,0,0],[0,0,0,1,0],[0,0,0,1,0],[1,0,0,0,0],[1,0,0,0,0],
[0,1,0,0,0],[1,0,0,0,0],[0,0,1,0,0],[0,0,1,0,0],[1,0,0,0,0],[0,0,0,1,0]])

n_ordered = res.shape[1]
coords = {'obs':np.arange(len(inp)), 'inp': np.arange(max(inp)+1), 'outp': np.arange(res.shape[1])}
with pm.Model(coords=coords) as mdl:

    mu = pm.Normal("intercept", sigma=3.5)[None]

    offset = pm.Normal("offset", dims=('inp'),
                 transform=pm.distributions.transforms.ZeroSumTransform([0]))
        
    scale = 3.5*pm.HalfStudentT("scale",nu=5)
    mu += (scale*offset)[inp]

    phi_delta = pm.Dirichlet(f'phi_diffs', [1.0]*(n_ordered-1) )
    phi = pt.concatenate([[0], pt.cumsum(phi_delta)])
    s_mu = pm.Normal(f"stereotype_intercept", size=n_ordered, 
                                transform=pm.distributions.transforms.ZeroSumTransform([-1]))
    fprobs = pm.math.softmax(s_mu[None,:] + phi[None,:]*mu[:,None], axis=-1)

    pm.Multinomial("y_res", p=fprobs, n=np.ones(len(inp)), observed=res, dims=('obs','outp'))

    # Seed 2 gives lbfgs_status, 4 gives "failed to iterate"
    idata = pme.fit_pathfinder(jitter=12,random_seed=np.random.default_rng(2))

[velochy]

Also, for the record, I have a larger model that goes into what seems like an infinite loop (with the default jitter value of 2)
Or at least I waited for 8 hours, by which NUTS would have completed the model, and it failed to perform a single path.
This is a rather large model though, so extracting it as a standalone code will be hard. But I'm thinking maybe pickling a model object will work, and allow you to load and run it?

[aphc14]

I had copy-pasted code from an example that set jitter to 12... which clearly seems to be too large of a value for my data. Reducing it down to the default made the models sample nicely.

Thanks so much @velochy ! My initial thoughts are that setting jitter to a very large value shouldn't be causing "BUG: Failed to iterate" because high jitters should trigger another Exception, which increases the count for INIT_FAILED when printing out the final results. The error could be due to large jitter + <something_else>. But I'd need to play around with your modelling scenario and understand how L-BFGS is behaving to figure out what that <something_else> is.

[aphc14]

Also, for the record, I have a larger model that goes into what seems like an infinite loop (with the default jitter value of 2) Or at least I waited for 8 hours, by which NUTS would have completed the model, and it failed to perform a single path. This is a rather large model though, so extracting it as a standalone code will be hard. But I'm thinking maybe pickling a model object will work, and allow you to load and run it?

Just to confirm, when you waited 8 hours for PF to finish, did you get any outputs? Or was it completely stuck at 0/num_paths completed in the progress bar? How long did NUTS take to finish?

I'm curious to know why PF is taking much longer than NUTS for your scenario. Would it be okay if you submit a separate issue with some extra details? Pickling the model + idata is an excellent idea. Found this link on how to pickle a pymc model: https://discourse.pymc.io/t/how-save-pymc-v5-models/13022/3

Not to worry if you cannot pickle the model + idata. I can just rerun your scenario on my end.

Thanks again @velochy !

[velochy]

Just to confirm, when you waited 8 hours for PF to finish, did you get any outputs? Or was it completely stuck at 0/num_paths completed in the progress bar? How long did NUTS take to finish?

Completely stuck. As for how long NUTS takes, I don't remember exactly. Started it again and after a few samples it shows ~5h, so probably in that ballpark. But in most cases, pathfinder is an order of magnitude faster so this looks very off.

Not to worry if you cannot pickle the model + idata. I can just rerun your scenario on my end.

Not so easy, unfortunately. We are often working with private data, and I am building a proprietary modeling stack on top of it, so I cannot really share neither the code that creates the model nor the data I use to run it, unless you are interested in signing a few NDAs :D

But if pickling works, I can make up some fake data similar to real data (or use some of the less private datasets) and pickle to model that comes out of the generator and that workflow should be relatively straightforward to file the reports for you.
I'll do a separate issue once I have some time to try doing it.

[aphc14]

If NUTS took ~5h, then something is very off with PF making zero progress after 8 hours. I'd be happy to investigate ^_^ But how would be tricky in this case.

After reading some horror stories about unpickling shared pickled files onto a personal machine, I think it'll be better to provide code snippets instead that provide a minimal reproducing example. It could be a significantly smaller version of your model that shows NUTS working and PF getting stuck or NUTS having superior speedup advantage over PF. Even one where NUTS takes ~15 minutes, and PF takes much longer or gets stuck. The data can be one that is randomly generated, as long as the model + data outputs show where PF is problematic.

Alternatively, you can try and add logging messages across the PF workflow but that could get quite messy though.

Not keen on any NDAs! 😂

[aphc14]

hi @velochy , appreciate you sharing this issue. I have been able to figure out the cause and the committed changes above ensures that the LBFGS iterations and PF algorithm is more robust to unstable trajectories 👍