hahnec / torchimize Goto Github PK
View Code? Open in Web Editor NEWnumerical optimization in pytorch
Home Page: https://hahnec.github.io/torchimize/
License: GNU General Public License v3.0
numerical optimization in pytorch
Home Page: https://hahnec.github.io/torchimize/
License: GNU General Public License v3.0
Comment out the the jacobian function input from the parallel fitter, e.g., here:
torchimize/tests/unit_test_parallel.py
Line 234 in 2ccbec4
and the test fails with:
======================================================================
ERROR: test_lma_emg_conditions (__main__.ParallelOptimizationTest)
----------------------------------------------------------------------
Traceback (most recent call last):
File "<ipython-input-59-7536f068804a>", line 215, in test_lma_emg_conditions
coeffs = lsq_lma_parallel(
File "/usr/local/lib/python3.10/dist-packages/torchimize/functions/parallel/lma_fun_parallel.py", line 95, in lsq_lma_parallel
p, f, g, H = newton_step_parallel(p, fun, jac_fun, wvec)
File "/usr/local/lib/python3.10/dist-packages/torchimize/functions/parallel/newton_parallel.py", line 40, in newton_step_parallel
gc = torch.einsum('bcnp,bcnp->bcp', j, f[..., None])
File "/usr/local/lib/python3.10/dist-packages/torch/functional.py", line 378, in einsum
return _VF.einsum(equation, operands) # type: ignore[attr-defined]
RuntimeError: einsum(): the number of subscripts in the equation (4) does not match the number of dimensions (5) for operand 0 and no ellipsis was given
Version 0.0.14 installed from pypi using pip
With chain rule one often sees many common sub-expressions in the calculation of the function and its derivatives. It would save time to compute the two of these together rather than having one call to the function and one call to the jacobian.
Fit results for any given fit depend on which other fits are happening in parallel.
One issue is the stopping condition, which applies to all fits simultaneously:
torchimize/torchimize/functions/parallel/lma_fun_parallel.py
Lines 108 to 111 in 2ccbec4
An individual fit may take more or fewer steps in batch than it would have done if it were performed alone. More because it continues to update even when it has converged because other fits in the block have not yet converged. Fewer because the badness in the residuals of any one fit is spread across the residuals in all fits.
A two-evaluation would help, using sum/max over the data dimensions to produce a vector of individual stopping criteria. You can then stop if all criteria are met. This will address the case of stopping too soon for some datasets.
For the other case you can suppress the update of the parameters for those fits that have already converged. This doesn't save any computation time, and will overall lead to worse fits (assuming more steps are always better when fitting), but it does mean that looping over single fits should in principle give the same results as performing fits in batch.
Hi,
thank you for your library, I'd like to use the LM optimizer in my PINN which is set up as a class based on pytorch. Could you share some examples of the code where you use the LM as an optimizer?
Thanks again!
Replace line with h = -l * torch.linalg.lstsq(j, f, rcond=None, driver=None)[0]
and skip the other matmuls and transposes.
Attempting to fit a set of simple gaussians I'm encountering frequent torch lstsq errors, presumably because the matrix is ill-conditioned near the minimum.
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
[<ipython-input-4-bde2f34b56b2>](https://localhost:8080/#) in <cell line: 69>()
67 extra = dict(meth = "mar", ftol = 0, ptol = 1e-8, gtol = 0)
68
---> 69 coeffs = opt(
70 p = p0,
71 function = lm_fun,
[/usr/local/lib/python3.10/dist-packages/torchimize/functions/parallel/lma_fun_parallel.py](https://localhost:8080/#) in lsq_lma_parallel(p, function, jac_function, args, wvec, ftol, ptol, gtol, tau, meth, rho1, rho2, beta, gama, max_iter)
96 D = lm_dg_step(H, D)
97 Hu = H+u[:, None, None]*D
---> 98 h = -torch.linalg.lstsq(Hu.double(), g.double(), rcond=None, driver=None)[0].to(dtype=p.dtype)
99 f_h = fun(p+h)
100 rho_nom = torch.einsum('bcp,bci->bc', f, f).sum(1) - torch.einsum('bcp,bci->bc', f_h, f_h).sum(1)
RuntimeError: false INTERNAL ASSERT FAILED at "../aten/src/ATen/native/BatchLinearAlgebra.cpp":1538, please report a bug to PyTorch. torch.linalg.lstsq: (Batch element 0): Argument 4 has illegal value. Most certainly there is a bug in the implementation calling the backend library.
torch version: 2.0.1+cu118
The following runs in a cell in colab. Maybe need to run several times to get it to fail, or increase n. Other optimizers fail similarly.
[Edit: Δy values were too small in some cases. Code modified to remove that problem.]
%pip install torchimize
n = 10
import torch
import numpy as np
sqrt2π = torch.sqrt(torch.tensor(2*torch.pi))
def lm_fx(p, x):
A, μ, σ = p.T
x, A, μ, σ = x[None, :], A[:, None], μ[:, None], σ[:, None]
fx = A*torch.exp(((x - μ)/σ)**2/-2)
return fx
def lm_fun(p, x, y, dy=1):
fx = lm_fx(p, x)
cost = ((fx - y)/dy)**2
return torch.unsqueeze(cost, 1) # Add an empty multi-cost dimension
def lm_jac(p, x, y, dy=1):
A, μ, σ = p.T
x, A, μ, σ = x[None, :], A[:, None], μ[:, None], σ[:, None]
shift = x - μ
scale = shift/σ
scalesq = scale**2
G = torch.exp(scalesq/-2)
fx = A * G
dfdA = G
dfdμ = fx * (scale/σ)
dfdσ = fx * (scalesq/σ)
Δ = (fx - y)/dy
#cost = Δ**2
partials = torch.stack([dfdA, dfdμ, dfdσ], dim=2)
#print(partials.shape, Δ.shape)
jac = (2*Δ[..., None])*partials
return torch.unsqueeze(jac, 1) # Add an empty multi-cost dimension
def gen(n, noise=0):
x = torch.linspace(0, 50, 30)
mid = (x[0] + x[-1])/2
width = (x[-1] - x[0])
μ, σ, A = torch.randn(n)*width/3 + mid, torch.rand(n)*width/5, 10**(1+torch.rand(n))
minwidth = 5
σ[σ < minwidth] = minwidth
#A = torch.ones_like(A)
pars = torch.stack([A, μ, σ]).T
y = lm_fx(pars, x)
if noise > 0:
dy = noise*y
dy[dy < noise] = noise
y += torch.randn(*y.shape)*dy
else:
dy = 1
return x, y, dy, pars
x, y, dy, target = gen(n, noise=0.05)
# Generate initial guess from the data
np_y, np_x = y.numpy(), x.numpy()
μ = np_x[np_y.argmax(axis=1)]
# Cheat for now and use target value rather than finding FWHM/2.35
σ = target[:, 2].numpy()
A = np_y.max(axis=1)
p0 = torch.tensor(np.vstack([A, μ, σ]).T)
# Bigger cheat: Use the target values as the intial guess
#p0 = target
from torchimize.functions.parallel.lma_fun_parallel import lsq_lma_parallel as opt
extra = dict(meth = "mar", ftol = 0, ptol = 1e-8, gtol = 0)
coeffs = opt(
p = p0,
function = lm_fun,
jac_function = lm_jac,
args = (x.to(dtype=p0.dtype), y.to(dtype=p0.dtype), dy),
max_iter = 99,
**extra,
)
print("num steps", len(coeffs))
best = coeffs[-1]
from matplotlib import pyplot as plt
%matplotlib inline
fx = lm_fx(best, x).detach()
for k in range(min(y.shape[0], 10)):
h, = plt.plot(x, y[k], '.')
plt.plot(x, fx[k], '-', color=h.get_color())
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.