Talk:Chebyshev iteration
 | This article is rated Start-class on Wikipedia's content assessment scale. It is of interest to the following WikiProjects: | ||||||||||
| |||||||||||
Code implementation
[edit]Hi,
It seems the Matlab code implementation doesn't match what the reference [1] proposes. If it is supposed to take improvements from [2], it should specify which algorithm it is using since [2] proposes 6 algorithm variants.
My proposal is to provide the Matlab code that matches [1] and then if required add another section on the improvements proposed by [2]. This should hopefully make the article much clearer. Varagrawal (talk) 19:00, 23 September 2022 (UTC)
- I do not think we should provide Matlab code that matches [1]. The reason is that it appears to have a bug. The code that used to be on the wiki page was translated from [1]. I've just independently tried to translate [1] into Matlab. I've also taken the reference code provided by [1] (available at https://netlib.org/templates/matlab/cheby.m). All three of these do not converge at the correct rate. The expected convergence rate can be found in Iterative Methods, Gutknecht, 2008, Theorem 1.6.2. In contrast, the code that is currently on the wiki page converges at the expected rate. I find it very surprising that the algorithm and code from such a popular reference as [1] appears to have a bug in it, but it looks that way.
- I agree that the code on the wiki page should be a direct and clear translation of some reference. At present, it's something I created by starting with the old code based on [1] and finding a minimal fix by tracing through the math in [2]. I don't think it's even an implementation of any particular algorithm from [2]. I don't think this sort of original development is the right thing to do on Wikipedia. But, we need to find a reference that converges at the correct rate, free from whatever bug is in [1]. I don't have such a reference. Perhaps one of the algorithms in [2] actually works. If we cannot find a reference with correct code, perhaps deleting the code altogether is the right thing to do? EssexEdwards (talk) 03:29, 5 October 2022 (UTC)
- After taking another look, a direct translation of algorithm 6 from [2] (https://people.math.ethz.ch/~mhg/pub/Cheby-02ParComp.pdf) converges at the expected rate. The algorithm is almost identical to the current code on the wiki page. One difference is a tiny algebra manipulation. The other is that the wiki page supports a preconditioner and the reference [2] does not. Also, all of the variable names are different, but that's to be expected.
- How about we provide the Matlab code the matches Algorithm 6 from [2] with the one-line addition to support preconditioning, and point out that one-line difference very clearly?
- For reference, here's my translation of Algorithm 6. It still needs some cleanup.
- EssexEdwards (talk) 14:49, 5 October 2022 (UTC)
::% From https://people.math.ethz.ch/~mhg/pub/Cheby-02ParComp.pdf ::% Algorithm 6. ::function [x, r_norm_history] = versionAlgorithm6(A, b, x0, iterNum, lMin, lMax) :: alpha = (lMax + lMin) / 2; :: c = (lMax - lMin) / 2; :: eta = -alpha/c; :: :: x = x0; :: r = b - A * x; :: v = 0*r; :: :: r_norm_history = norm(r) * ones(iterNum+1,1); :: for i = 1:iterNum :: if (i == 1) :: phi = 0; :: omega = 1/alpha; :: elseif (i == 2) :: phi = -0.5*(c/alpha)^2; :: omega = 1.0/(alpha - c^2/(2*alpha)); :: else :: phi = -(c/2)^2 * omega^2; :: omega = 1.0/(alpha - (c/2)^2*omega); :: end :: v = r - v * phi; :: x = x + v * omega; :: :: r = b - A * x; :: r_norm_history(i+1) = norm(r); :: :: if (norm(r) < 1e-15), break; end; % stop if necessary :: end ::end ::