Talk:Chebyshev iteration

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 [https://people.math.ethz.ch/~mhg/unt/SWNLA/itmethSWNLA08.pdf 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.

::

::% 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

:: EssexEdwards (talk) 14:49, 5 October 2022 (UTC)