celerite2.backprop.LinAlgError: failed to factorize or solve matrix

**Describe the bug**
pm.Model() returns

> celerite2.backprop.LinAlgError: failed to factorize or solve matrix

**To Reproduce**
While following the [Fitting TESS data](https://gallery.exoplanet.codes/tutorials/tess/) tutorial, I try to fit LHS 3844. I used my own light curve (30 min TESS FFI), and the only things I changed in the code are the stellar mass and radius prior:
```
# Stellar parameters from Parallax +Benedict et al. (2016)
M_star= 0.151, 0.014
# Parallax +Mann et al. (2015)
R_star = 0.189, 0.006
```
The period is 0.4626. I had to change the BLS range to include this period. But BLS worked fine. I suspect this short period is a problem.

**Your setup (please complete the following information):**
- Version of exoplanet: 0.5.2
- Operating system: WLS Ubuntu 20.04
- Python version & installation method (pip, conda, etc.): 3.8.5, o

**Additional context**
Full error message:

```
optimizing logp for variables: [log_sigma_gp, log_rho_gp, log_sigma_lc, ecs, log_depth, b, log_period, t0, r_star, m_star, u_star, mean]

array: [ 12.90461679 28.42071018 0.17878725 0.31305992 0.53748097
-6.02102727 9.10216664 0.21911012 -23.82860354 -2.70012843
-4.1448473 0.27730417 0.13437652 35.9584247 ]
point: {'mean': array(35.9584247), 'u_star_quadlimbdark__': array([0.27730417, 0.13437652]), 'm_star_interval__': array(-4.1448473), 'r_star_interval__': array(-2.70012843), 't0': array(-23.82860354), 'log_period': array(0.21911012), 'b_interval__': array(9.10216664), 'log_depth': array(-6.02102727), 'ecs_unitdisk+interval__': array([0.31305992, 0.53748097]), 'log_sigma_lc': array(0.17878725), 'log_rho_gp': array(28.42071018), 'log_sigma_gp': array(12.90461679)}
Traceback (most recent call last):
File "/home/tehan/.local/lib/python3.8/site-packages/theano/compile/function/types.py", line 974, in __call__
self.fn()
File "/home/tehan/.local/lib/python3.8/site-packages/theano/graph/op.py", line 476, in rval
r = p(n, [x[0] for x in i], o)
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py", line 123, in perform
func(*args)
celerite2.backprop.LinAlgError: failed to factorize or solve matrix

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/home/tehan/.local/lib/python3.8/site-packages/pymc3_ext/optim.py", line 213, in __call__
res = self.func(
File "/home/tehan/.local/lib/python3.8/site-packages/theano/compile/function/types.py", line 987, in __call__
raise_with_op(
File "/home/tehan/.local/lib/python3.8/site-packages/theano/link/utils.py", line 508, in raise_with_op
raise exc_value.with_traceback(exc_trace)
File "/home/tehan/.local/lib/python3.8/site-packages/theano/compile/function/types.py", line 974, in __call__
self.fn()
File "/home/tehan/.local/lib/python3.8/site-packages/theano/graph/op.py", line 476, in rval
r = p(n, [x[0] for x in i], o)
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py", line 123, in perform
func(*args)
celerite2.backprop.LinAlgError: failed to factorize or solve matrix
Apply node that caused the error: _CeleriteOp{name='factor_fwd', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0)
Toposort index: 135
Inputs types: [TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, matrix), TensorType(float64, matrix)]
Inputs shapes: [(1821,), (2,), (1821,), (1821, 2), (1821, 2)]
Inputs strides: [(8,), (8,), (8,), (16, 8), (16, 8)]
Inputs values: ['not shown', array([2.01700636e-12, 2.01700636e-12]), 'not shown', 'not shown', 'not shown']
Outputs clients: [[_CeleriteOp{name='factor_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, _CeleriteOp{name='factor_fwd', quiet=False}.2, Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}.0, _CeleriteOp{name='solve_lower_rev', quiet=False}.3), Elemwise{true_div,no_inplace}(InplaceDimShuffle{0}.0, _CeleriteOp{name='factor_fwd', quiet=False}.0), Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}(TensorConstant{(1,) of 0.5}, _CeleriteOp{name='factor_fwd', quiet=False}.0, TensorConstant{(1,) of 0.5}, TensorConstant{(1,) of -1.0}, Elemwise{Sqr}[(0, 0)].0), Elemwise{TrueDiv}[(0, 0)](Elemwise{Sqr}[(0, 0)].0, _CeleriteOp{name='factor_fwd', quiet=False}.0), Elemwise{Log}[(0, 0)](_CeleriteOp{name='factor_fwd', quiet=False}.0)], [_CeleriteOp{name='solve_lower_fwd', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, Elemwise{Composite{(i0 - ((i1 * i2) + i3))}}[(0, 2)].0), _CeleriteOp{name='solve_lower_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, Elemwise{Composite{(i0 - ((i1 * i2) + i3))}}[(0, 2)].0, _CeleriteOp{name='solve_lower_fwd', quiet=False}.0, _CeleriteOp{name='solve_lower_fwd', quiet=False}.1, IncSubtensor{InplaceInc;::, int64}.0), _CeleriteOp{name='factor_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, _CeleriteOp{name='factor_fwd', quiet=False}.2, Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}.0, _CeleriteOp{name='solve_lower_rev', quiet=False}.3)], [_CeleriteOp{name='factor_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, _CeleriteOp{name='factor_fwd', quiet=False}.2, Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}.0, _CeleriteOp{name='solve_lower_rev', quiet=False}.3)]]

Backtrace when the node is created(use Theano flag traceback__limit=N to make it longer):
File "", line 141, in
model0, map_soln0, extras0 = build_model()
File "", line 98, in build_model
gp = GaussianProcess(kernel, t=x[mask], yerr=tt.exp(log_sigma_lc))
File "/usr/local/lib/python3.8/dist-packages/celerite2/core.py", line 210, in __init__
self.compute(t, **kwargs)
File "/usr/local/lib/python3.8/dist-packages/celerite2/core.py", line 316, in compute
self._do_compute(quiet)
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/celerite2.py", line 83, in _do_compute
self._d, self._W, _ = ops.factor(
File "/home/tehan/.local/lib/python3.8/site-packages/theano/graph/op.py", line 250, in __call__
node = self.make_node(*inputs, **kwargs)
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py", line 83, in make_node
otypes = [
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py", line 84, in
tt.TensorType(

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
File ~/.local/lib/python3.8/site-packages/theano/compile/function/types.py:974, in Function.__call__(self, *args, **kwargs)
972 try:
973 outputs = (
--> 974 self.fn()
975 if output_subset is None
976 else self.fn(output_subset=output_subset)
977 )
978 except Exception:

File ~/.local/lib/python3.8/site-packages/theano/graph/op.py:476, in Op.make_py_thunk..rval(p, i, o, n)
475 def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 476 r = p(n, [x[0] for x in i], o)
477 for o in node.outputs:

File /usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py:123, in _CeleriteOp.perform(self, node, inputs, outputs)
122 try:
--> 123 func(*args)
124 except backprop.LinAlgError:

LinAlgError: failed to factorize or solve matrix

During handling of the above exception, another exception occurred:

LinAlgError Traceback (most recent call last)
Input In [6], in
131 extras = dict(
132 zip(
133 ["light_curves", "gp_pred"],
134 pmx.eval_in_model([light_curves, gp.predict(resid)], map_soln),
135 )
136 )
138 return model, map_soln, extras
--> 141 model0, map_soln0, extras0 = build_model()

Input In [6], in build_model(mask, start)
125 map_soln = pmx.optimize(start=map_soln, vars=[mean])
126 map_soln = pmx.optimize(
127 start=map_soln, vars=[log_sigma_lc, log_sigma_gp, log_rho_gp]
128 )
--> 129 map_soln = pmx.optimize(start=map_soln)
131 extras = dict(
132 zip(
133 ["light_curves", "gp_pred"],
134 pmx.eval_in_model([light_curves, gp.predict(resid)], map_soln),
135 )
136 )
138 return model, map_soln, extras

File ~/.local/lib/python3.8/site-packages/pymc3_ext/optim.py:119, in optimize(start, vars, return_info, verbose, progress, maxeval, model, **kwargs)
116 kwargs["jac"] = True
118 try:
--> 119 info = minimize(objective, x0, **kwargs)
120 except (KeyboardInterrupt, StopIteration):
121 info = None

File ~/.local/lib/python3.8/site-packages/scipy/optimize/_minimize.py:618, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
616 return _minimize_cg(fun, x0, args, jac, callback, **options)
617 elif meth == 'bfgs':
--> 618 return _minimize_bfgs(fun, x0, args, jac, callback, **options)
619 elif meth == 'newton-cg':
620 return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
621 **options)

File ~/.local/lib/python3.8/site-packages/scipy/optimize/optimize.py:1235, in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, finite_diff_rel_step, **unknown_options)
1232 pk = -np.dot(Hk, gfk)
1233 try:
1234 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
-> 1235 _line_search_wolfe12(f, myfprime, xk, pk, gfk,
1236 old_fval, old_old_fval, amin=1e-100, amax=1e100)
1237 except _LineSearchError:
1238 # Line search failed to find a better solution.
1239 warnflag = 2

File ~/.local/lib/python3.8/site-packages/scipy/optimize/optimize.py:1005, in _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval, **kwargs)
991 """
992 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
993 suitable step length is not found, and raise an exception if a
(...)
1000
1001 """
1003 extra_condition = kwargs.pop('extra_condition', None)
-> 1005 ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
1006 old_fval, old_old_fval,
1007 **kwargs)
1009 if ret[0] is not None and extra_condition is not None:
1010 xp1 = xk + ret[0] * pk

File ~/.local/lib/python3.8/site-packages/scipy/optimize/linesearch.py:96, in line_search_wolfe1(f, fprime, xk, pk, gfk, old_fval, old_old_fval, args, c1, c2, amax, amin, xtol)
92 return np.dot(gval[0], pk)
94 derphi0 = np.dot(gfk, pk)
---> 96 stp, fval, old_fval = scalar_search_wolfe1(
97 phi, derphi, old_fval, old_old_fval, derphi0,
98 c1=c1, c2=c2, amax=amax, amin=amin, xtol=xtol)
100 return stp, fc[0], gc[0], fval, old_fval, gval[0]

File ~/.local/lib/python3.8/site-packages/scipy/optimize/linesearch.py:172, in scalar_search_wolfe1(phi, derphi, phi0, old_phi0, derphi0, c1, c2, amax, amin, xtol)
170 if task[:2] == b'FG':
171 alpha1 = stp
--> 172 phi1 = phi(stp)
173 derphi1 = derphi(stp)
174 else:

File ~/.local/lib/python3.8/site-packages/scipy/optimize/linesearch.py:84, in line_search_wolfe1..phi(s)
82 def phi(s):
83 fc[0] += 1
---> 84 return f(xk + s*pk, *args)

File ~/.local/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py:249, in ScalarFunction.fun(self, x)
247 if not np.array_equal(x, self.x):
248 self._update_x_impl(x)
--> 249 self._update_fun()
250 return self.f

File ~/.local/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py:233, in ScalarFunction._update_fun(self)
231 def _update_fun(self):
232 if not self.f_updated:
--> 233 self._update_fun_impl()
234 self.f_updated = True

File ~/.local/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py:137, in ScalarFunction.__init__..update_fun()
136 def update_fun():
--> 137 self.f = fun_wrapped(self.x)

File ~/.local/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py:134, in ScalarFunction.__init__..fun_wrapped(x)
130 self.nfev += 1
131 # Send a copy because the user may overwrite it.
132 # Overwriting results in undefined behaviour because
133 # fun(self.x) will change self.x, with the two no longer linked.
--> 134 return fun(np.copy(x), *args)

File ~/.local/lib/python3.8/site-packages/scipy/optimize/optimize.py:74, in MemoizeJac.__call__(self, x, *args)
72 def __call__(self, x, *args):
73 """ returns the the function value """
---> 74 self._compute_if_needed(x, *args)
75 return self._value

File ~/.local/lib/python3.8/site-packages/scipy/optimize/optimize.py:68, in MemoizeJac._compute_if_needed(self, x, *args)
66 if not np.all(x == self.x) or self._value is None or self.jac is None:
67 self.x = np.asarray(x).copy()
---> 68 fg = self.fun(x, *args)
69 self.jac = fg[1]
70 self._value = fg[0]

File ~/.local/lib/python3.8/site-packages/pymc3_ext/optim.py:103, in optimize..objective(vec)
101 nonlocal neval
102 neval += 1
--> 103 nll, grad = wrapper(vec)
104 if progressbar:
105 progressbar.comment = "logp = {0:.3e}".format(-nll)

File ~/.local/lib/python3.8/site-packages/pymc3_ext/optim.py:213, in ModelWrapper.__call__(self, vec)
211 def __call__(self, vec):
212 try:
--> 213 res = self.func(
214 *get_args_for_theano_function(
215 self.bij.rmap(vec), model=self.model
216 )
217 )
218 except Exception:
219 import traceback

File ~/.local/lib/python3.8/site-packages/theano/compile/function/types.py:987, in Function.__call__(self, *args, **kwargs)
985 if hasattr(self.fn, "thunks"):
986 thunk = self.fn.thunks[self.fn.position_of_error]
--> 987 raise_with_op(
988 self.maker.fgraph,
989 node=self.fn.nodes[self.fn.position_of_error],
990 thunk=thunk,
991 storage_map=getattr(self.fn, "storage_map", None),
992 )
993 else:
994 # old-style linkers raise their own exceptions
995 raise

File ~/.local/lib/python3.8/site-packages/theano/link/utils.py:508, in raise_with_op(fgraph, node, thunk, exc_info, storage_map)
505 warnings.warn(f"{exc_type} error does not allow us to add extra error message")
506 # Some exception need extra parameter in inputs. So forget the
507 # extra long error message in that case.
--> 508 raise exc_value.with_traceback(exc_trace)

File ~/.local/lib/python3.8/site-packages/theano/compile/function/types.py:974, in Function.__call__(self, *args, **kwargs)
971 t0_fn = time.time()
972 try:
973 outputs = (
--> 974 self.fn()
975 if output_subset is None
976 else self.fn(output_subset=output_subset)
977 )
978 except Exception:
979 restore_defaults()

File ~/.local/lib/python3.8/site-packages/theano/graph/op.py:476, in Op.make_py_thunk..rval(p, i, o, n)
475 def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 476 r = p(n, [x[0] for x in i], o)
477 for o in node.outputs:
478 compute_map[o][0] = True

File /usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py:123, in _CeleriteOp.perform(self, node, inputs, outputs)
120 func = getattr(driver, self.name)
122 try:
--> 123 func(*args)
124 except backprop.LinAlgError:
125 if not self.quiet:

LinAlgError: failed to factorize or solve matrix
Apply node that caused the error: _CeleriteOp{name='factor_fwd', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0)
Toposort index: 135
Inputs types: [TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, matrix), TensorType(float64, matrix)]
Inputs shapes: [(1821,), (2,), (1821,), (1821, 2), (1821, 2)]
Inputs strides: [(8,), (8,), (8,), (16, 8), (16, 8)]
Inputs values: ['not shown', array([2.01700636e-12, 2.01700636e-12]), 'not shown', 'not shown', 'not shown']
Outputs clients: [[_CeleriteOp{name='factor_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, _CeleriteOp{name='factor_fwd', quiet=False}.2, Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}.0, _CeleriteOp{name='solve_lower_rev', quiet=False}.3), Elemwise{true_div,no_inplace}(InplaceDimShuffle{0}.0, _CeleriteOp{name='factor_fwd', quiet=False}.0), Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}(TensorConstant{(1,) of 0.5}, _CeleriteOp{name='factor_fwd', quiet=False}.0, TensorConstant{(1,) of 0.5}, TensorConstant{(1,) of -1.0}, Elemwise{Sqr}[(0, 0)].0), Elemwise{TrueDiv}[(0, 0)](Elemwise{Sqr}[(0, 0)].0, _CeleriteOp{name='factor_fwd', quiet=False}.0), Elemwise{Log}[(0, 0)](_CeleriteOp{name='factor_fwd', quiet=False}.0)], [_CeleriteOp{name='solve_lower_fwd', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, Elemwise{Composite{(i0 - ((i1 * i2) + i3))}}[(0, 2)].0), _CeleriteOp{name='solve_lower_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, Elemwise{Composite{(i0 - ((i1 * i2) + i3))}}[(0, 2)].0, _CeleriteOp{name='solve_lower_fwd', quiet=False}.0, _CeleriteOp{name='solve_lower_fwd', quiet=False}.1, IncSubtensor{InplaceInc;::, int64}.0), _CeleriteOp{name='factor_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, _CeleriteOp{name='factor_fwd', quiet=False}.2, Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}.0, _CeleriteOp{name='solve_lower_rev', quiet=False}.3)], [_CeleriteOp{name='factor_rev', quiet=False}(TensorConstant{[-24.07638...07638784]}, Join.0, Alloc.0, Join.0, Join.0, _CeleriteOp{name='factor_fwd', quiet=False}.0, _CeleriteOp{name='factor_fwd', quiet=False}.1, _CeleriteOp{name='factor_fwd', quiet=False}.2, Elemwise{Composite{((i0 / i1) + ((i2 * i3 * i4) / sqr(i1)))}}.0, _CeleriteOp{name='solve_lower_rev', quiet=False}.3)]]

Backtrace when the node is created(use Theano flag traceback__limit=N to make it longer):
File "", line 141, in
model0, map_soln0, extras0 = build_model()
File "", line 98, in build_model
gp = GaussianProcess(kernel, t=x[mask], yerr=tt.exp(log_sigma_lc))
File "/usr/local/lib/python3.8/dist-packages/celerite2/core.py", line 210, in __init__
self.compute(t, **kwargs)
File "/usr/local/lib/python3.8/dist-packages/celerite2/core.py", line 316, in compute
self._do_compute(quiet)
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/celerite2.py", line 83, in _do_compute
self._d, self._W, _ = ops.factor(
File "/home/tehan/.local/lib/python3.8/site-packages/theano/graph/op.py", line 250, in __call__
node = self.make_node(*inputs, **kwargs)
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py", line 83, in make_node
otypes = [
File "/usr/local/lib/python3.8/dist-packages/celerite2/theano/ops.py", line 84, in
tt.TensorType(

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.
```

Thank you for your help!