The coolest descent to find a solution to a linear system with a Hilbert matrix

Question

I use the steepest descent method to figure out a solution to a linear system with a 5x5 Hilbert matrix. I believe the code is good in the sense that it gives me the correct answer.

My problem is that:

• I think too many iterations are suitable for the correct answer. I believe that maybe something was missing in the algorithm, but I'm not sure at the moment.

• I’m not sure if this is the most efficient way to implement the algorithm, and besides, it’s a bit confusing on what tol should be chosen.

Any understanding of these issues will be appreciated (especially 1.). Thanks!

% Method of Steepest Descent with tol 10^-6 h = hilb(5); %Hilbert 5x5 matrix b = [1;1;1;1;1]; %solution matrix solution = zeros(d,1); %Initialization residual = h*solution - b; tol = 10^(-6) count = 0; while residual'*residual > tol; roe = (residual'*residual)/(residual'*h*residual); solution = solution - roe*residual; residual = h*solution - b; count = count + 1; end count solution %Method of Steepest Descent with tol 10^-12 solution = zeros(d,1); residual = h*solution - b; tol = 10^(-12) count = 0; while residual'*residual > tol; roe = (residual'*residual)/(residual'*h*residual); solution = solution - roe*residual; residual = residual - roe*h*residual; count = count + 1; end count solution %another_solution = invhilb(5)*b %Check for solution 

Answer 1

I have no knowledge to solve your problem from a mathematical point of view, but from a programming point of view there is a point that I could notice.

Indeed, you are right. It takes too many iterations until it reaches the result:

 Elapsed time is 6.487824 seconds. count = 292945 

When you determine the step size to achieve the final result, the step length must be optimized. Otherwise, you either approach the answer too slowly, or go through the optimal answer several times and zigzag around it, because the length of your step is too long.

To optimize the step size, I first create a function according to your script (plus some minor modifications):

 function [solution, count, T] = SDhilb(d, step) h = hilb(d); tic b = ones(d, 1); solution = zeros(d, 1); residual = h * solution - b; res2 = residual' * residual; tol = 10^(-6); count = 0; while res2 > tol; roe = res2/(residual' * h * residual); solution = solution - step * roe * residual; % here the step size appears residual = h * solution - b; count = count + 1; res2 = residual' * residual; % let calculate this once per iteration end T = toc; 

Now, using this function for a range of steps (0.1: 0.1: 1.3) and a pair of Hilbert matrices (d = 2, 3, 4, 5), it is obvious that 1 not an effective step size:

Instead, step = 0.9 seems a lot more efficient. let's see how effective it is in the case of hilb(5) :

 [result, count, T] = SDhilb(5, 0.9) result = 3.1894 -85.7689 481.4906 -894.8742 519.5830 count = 1633 T = 0.0204 

This is more than two orders of magnitude faster.

Similarly, you can try different tol values ​​to see how dramatic this can affect speed. In this case, there is no optimal value: the smaller the tolerance, the more time you need to wait.

Answer 2

It looks like you implemented the algorithm correctly (steep descent / gradient descent with exact linear search to minimize a convex quadratic function).

Convergence is slow because the problem is poorly conditioned: the Hilbert matrix that you are considering has a number of states above 400,000. Gradient descent is known to be slow when the problem is poorly conditioned.

Given instead a problem with good convention, for example, by adding identity to the Hilbert matrix (h = hilb (5) + eye (5)), the same code ends after 7 iterations (and the condition number for this matrix is ​​less than 3).