Hello,
I'm a co-developer of the package chemotools and for this I wanted to implement a deconvolution of spectra blurred by a lineshape function. I stumbled over your publication Gaussian Process Deconvolution. Since I'm not a native mathematician, I wanted to refer to the code provided for the publication to gain some further understanding.
When looking at the code, I've found some possible improvements that could be made:
- type hinting: Python allows for type hints to indicate the user which type is expected by a function/method. While these types are not enforced (unlike typed languages such as C, Fortran, or Rust), it makes it way easier to work with the code
- comments and docstrings: the code has almost no comments on why something is done, e.g.,
Gram = Spec_Mix(self.times_y,self.times_y,self.gamma,self.theta,self.sigma) + 1e-8*np.eye(self.Ny)
1e-8*np.eye(self.Ny) increases the numerical stability for the following solve could be really helpful. Docstrings for classes, functions, and methods point in a similar direction.
- adhering to Python conventions and using descriptive names: some functions are written with capital letters, e.g.,
Spec_Mix, SE, which is not consistent with the Python conventions. Something similar applies to the class names because they should contain capital letters. On top of that, the naming is often not really descriptive. What does SE do? It computes the squared exponential kernel, then why not name it something like squared_exponential?
- removing dead code: a lot of code is commented out and not used anymore.
- being explicit: the line
Gram = self.kernel_sum(self.times_y, self.times_y, params[:-1], weights, self.weight_step)
gram = self.kernel_sum(
x=self.times_y,
y=self.times_y,
params=params[:-1],
weights=weights,
weight_step=self.weight_step,
)self.kernel_sum and especially the order of its input arguments by heart. Now, everything is clear and due to the assignment to keywords, no order can be confused and lead to wrong output that might slip through unnoticed if the wrong order does not lead to a crash.
- formatting and linting: tools like
black and ruff help in keeping the code base clean and also for some of the previous points. Setting them up is almost no effort.
Sticking to 1 - 6 would be the first step towards a state that allows for actually improving the functionality of the code:
-
many critical parts are painfully slow, even though there are no complex computations: the reason is that they are wasting a significant amount of memory and runtime, e.g.,
(sign, logdet) = np.linalg.slogdet(K)
return 0.5*( Y.T@np.linalg.solve(K,Y) + logdet + self.Ny*np.log(2*np.pi))
where the from-scratch-inversion of K follows the computation of its log-determinant. This could be combined by using a Cholesky factorization or Pivoted LU-decomposition of K and then computing both the log-determinant and the subsequent inversion from it in virtually one go.
Another example is
dlogp_dsigma = sigma * 0.5*np.trace(H@dKdsigma)
A matrix product does never need to be fully formed if only the trace, and thus the main diagonal, has to be computed. Selectively evaluating the matrix product reduces the code to
dlogp_dsigma = sigma * 0.5*(H * dKdsigma.transpose()).sum()
which only uses a tiny (almost vanishing) amount of compute compared to the full product.
Another line offering room for improvement is:
Gram = Spec_Mix(self.times_y,self.times_y,self.gamma,self.theta,self.sigma) + 1e-8*np.eye(self.Ny)
where the main diagonal is updated by adding a matrix of mostly zeros. For Ny=100, only 1% of the added elements will be non-zero and this reduces to 0.1% already for Ny=1_000. On top of that, the overhead to initialise a diagonal matrix - where the main diagonal elements are quite badly scattered over the memory by nature - is significant.
One final example is:
def outersum(a,b):
return np.outer(a,np.ones_like(b))+np.outer(np.ones_like(a),b)for which NumPy already offers some powerful alternative. This line of code multiplies a vector with a lot of ones to reshape it. This is highly inefficient because the outer sum can also be achieved by using NumPy's broadcasting:
def outer_sum(a, b):
return a[::, np.newaxis] + b[np.newaxis, ::]Making custom logics for problems that so many others must have had already is both inefficient and also error-prone.
-
some parts are inconsistent: some methods behave the same, but yet with subtle differences, e.g.,
def neg_log_likelihood(self):
Y = self.y
Gram = Spec_Mix(self.times_y,self.times_y,self.gamma,self.theta,self.sigma) + 1e-8*np.eye(self.Ny)
K = Gram + self.sigma_n**2*np.eye(self.Ny)
(sign, logdet) = np.linalg.slogdet(K)
return 0.5*( Y.T@np.linalg.solve(K,Y) + logdet + self.Ny*np.log(2*np.pi))
def nlogp(self, hypers):
sigma = np.exp(hypers[0])
gamma = np.exp(hypers[1])
theta = np.exp(hypers[2])*0 #multiply by zero to force SE kernel (and not spectral mixture)
sigma_n = np.exp(hypers[3])
Y = self.y
Gram = Spec_Mix(self.times_y,self.times_y,gamma,theta,sigma)
K = Gram + sigma_n**2*np.eye(self.Ny) + 1e-5*np.eye(self.Ny)
(sign, logdet) = np.linalg.slogdet(K)
return 0.5*( Y.T@np.linalg.solve(K,Y) + logdet + self.Ny*np.log(2*np.pi))Can someone spot the difference? In neg_log_likelihood, 1e-8 is added to the main diagonal while in nlogp it's even 1e-5. Why is this like this? Wouldn't a class attribute like jitter (as this small summand is sometimes called) help keeping everything identical? Then, the user could even specify it and not rely on a hard-coded 1e-5 or 1e-8 to be a good choice.
This was what I've found from skimming the code, but I don't only want to complain. In contrary, I want to offer my help because in my opinion the code from the examples also massively questions the credibility of the publication, because all these facts make the code hardly maintainable/readable and also difficult to debug (independent of whether the code was meant to be on production quality level or not).
As an added bonus, the improvements will increase the audience and allow people to actually make use ot the publication code.
Please let me know if you are interested in a collaboration. Otherwise, I will work in a fork of the repository.
Kind regards
Hello,
I'm a co-developer of the package
chemotoolsand for this I wanted to implement a deconvolution of spectra blurred by a lineshape function. I stumbled over your publication Gaussian Process Deconvolution. Since I'm not a native mathematician, I wanted to refer to the code provided for the publication to gain some further understanding.When looking at the code, I've found some possible improvements that could be made:
1e-8*np.eye(self.Ny)increases the numerical stability for the followingsolvecould be really helpful. Docstrings for classes, functions, and methods point in a similar direction.Spec_Mix,SE, which is not consistent with the Python conventions. Something similar applies to the class names because they should contain capital letters. On top of that, the naming is often not really descriptive. What doesSEdo? It computes the squared exponential kernel, then why not name it something likesquared_exponential?self.kernel_sumand especially the order of its input arguments by heart. Now, everything is clear and due to the assignment to keywords, no order can be confused and lead to wrong output that might slip through unnoticed if the wrong order does not lead to a crash.blackandruffhelp in keeping the code base clean and also for some of the previous points. Setting them up is almost no effort.Sticking to 1 - 6 would be the first step towards a state that allows for actually improving the functionality of the code:
many critical parts are painfully slow, even though there are no complex computations: the reason is that they are wasting a significant amount of memory and runtime, e.g.,
where the from-scratch-inversion of
Kfollows the computation of its log-determinant. This could be combined by using a Cholesky factorization or Pivoted LU-decomposition ofKand then computing both the log-determinant and the subsequent inversion from it in virtually one go.Another example is
A matrix product does never need to be fully formed if only the trace, and thus the main diagonal, has to be computed. Selectively evaluating the matrix product reduces the code to
which only uses a tiny (almost vanishing) amount of compute compared to the full product.
Another line offering room for improvement is:
where the main diagonal is updated by adding a matrix of mostly zeros. For
Ny=100, only 1% of the added elements will be non-zero and this reduces to 0.1% already forNy=1_000. On top of that, the overhead to initialise a diagonal matrix - where the main diagonal elements are quite badly scattered over the memory by nature - is significant.One final example is:
for which NumPy already offers some powerful alternative. This line of code multiplies a vector with a lot of ones to reshape it. This is highly inefficient because the outer sum can also be achieved by using NumPy's broadcasting:
Making custom logics for problems that so many others must have had already is both inefficient and also error-prone.
some parts are inconsistent: some methods behave the same, but yet with subtle differences, e.g.,
Can someone spot the difference? In
neg_log_likelihood,1e-8is added to the main diagonal while innlogpit's even1e-5. Why is this like this? Wouldn't a class attribute likejitter(as this small summand is sometimes called) help keeping everything identical? Then, the user could even specify it and not rely on a hard-coded1e-5or1e-8to be a good choice.This was what I've found from skimming the code, but I don't only want to complain. In contrary, I want to offer my help because in my opinion the code from the examples also massively questions the credibility of the publication, because all these facts make the code hardly maintainable/readable and also difficult to debug (independent of whether the code was meant to be on production quality level or not).
As an added bonus, the improvements will increase the audience and allow people to actually make use ot the publication code.
Please let me know if you are interested in a collaboration. Otherwise, I will work in a fork of the repository.
Kind regards