Gauss-Newton Solver: Improper assignment with rectangular empty matrix

Question

I don't understand why I am getting this error message. I have implemented a Gauss-Newton solver to solve a system of linear equations.

It keeps saying improper assignment with rectangular empty matrix at the line "for i=1:m"

function [x, l] = GS(A, b, x0, TOL)
[m n] = size(A);
l = 1;
x = [0;0;0];
while (true)
    for i=1:m
        sum1 = 0; sum2 = 0;
    for j=1:i-1
        sum1 = sum1 + A(i,j)*x(j);
        for j=i+1:n
            sum2 = sum2 + A(i,j)*x(j);
        end
    end 
    x(i) = (-sum1-sum2+b(i))./A(i,j);    
end
if abs(norm(x) - norm(x0)) < TOL   
    break
end
x0 = x;
l = l + 1;
end

I have written an answer. I have also solved your other issue with your Jacobi solver. If you don't mind, take a look at my answer, and accept it if I have helped you... this one too! — rayryeng, Feb 11 '15 at 19:32
Thank you! I didn't know how to accept an answer but I see it now. Thank you for help on both of my codes! — user3681755, Feb 11 '15 at 19:46
No problem :) I actually made a slight mistake with my code. I swapped two variables within the two `for` loops within that bigger `for` loop using `i`. Take a look! It doesn't change the answer though... but it converges faster. — rayryeng, Feb 11 '15 at 19:48
Hmm, it seems that my Gauss Seidel code is computing the same as the Jacobi. It actually takes more iterations that Jacobi, which shouldn't happen because GS is an improved Jacobi method, correct? — user3681755, Feb 11 '15 at 20:37
That totally depends on the system setup. Sometimes, GS converges **faster** than Jacobi. Also, I made a small error when I wrote my initial answer to you. I had to swap two variables. — rayryeng, Feb 11 '15 at 20:39
Also, look here: http://stackoverflow.com/questions/24730993/jacobi-iteration-doesnt-end/24733027#24733027 . You may get a case where Jacobi **doesn't** converge while Gauss-Seidel does. It's the way you are conditioning the solution at each iteration. The subtle difference between using the current solution mixed with the previous solution in each iteration and only the previous solution can make a huge difference. — rayryeng, Feb 11 '15 at 20:40
If I could just get your help with one more thing... I am now using the SOR method and need to find the optimal weight factor. I think a good way to go about this is to run the SOR code with a number of omegas from 0 to 2, then store the number of iterations for each of these. Then I can see which iteration is the lowest and which omega it corresponds to. Being a novice programer, however, I am unsure how to go about this. — user3681755, Feb 11 '15 at 21:39
Make a new post. However, that's pretty simple. You'd simply loop over all of the omegas from 0 to 2, record how many iterations it took to solve, as well as the solution vector for that omega. Whichever gave you the minimum number of iterations is what the final solution vector would be. — rayryeng, Feb 11 '15 at 21:43

rayryeng · Accepted Answer · 2015-02-11T19:48:24.047

Three things wrong with your code.

Your for loop statements are mis-aligned. Specifically, you forgot to place an end keyword at the end of the first for loop within the i loop (i.e. for j = 1 : i-1)
You need to compute the values for Gauss-Newton using the current solution from j = 1, 2 up to i-1. The first for loop needs to use x, and the second for loop from j = i+1 up to n needs to use x0.
You need to divide by the correct diagonal coefficient when you solve for each value of x.

As such:

function [x, l] = GS(A, b, x0, TOL)
[m n] = size(A);
l = 1;
x = [0;0;0];
while (true)
    for i=1:m
        sum1 = 0; sum2 = 0;
        for j=1:i-1 %// CHANGE
            sum1 = sum1 + A(i,j)*x(j); %// CHANGE
        end
        for j=i+1:n %// CHANGE
            sum2 = sum2 + A(i,j)*x0(j);
        end
        x(i) = (-sum1-sum2+b(i))./A(i,i); %// CHANGE   
    end

    if abs(norm(x) - norm(x0)) < TOL   
        break
    end
x0 = x;
l = l + 1;
end

Example use:

format long;
A = [6 1 1; 1 5 3; 0 2 4]
b = [1 2 3].';
[x,i] = GS(A, b, [0;0;0], 1e-10)

x =

   0.048780487792648
  -0.085365853612062
   0.792682926806031


i =

    21

This means it took 21 iterations to achieve a solution with tolerance 1e-10. Compare this with MATLAB's built-in inverse:

x2 = A \ b

x2 =

   0.048780487804878
  -0.085365853658537
   0.792682926829268

As you can see, I specified a tolerance of 1e-10, which means we are guaranteed to have 10 decimal places of accuracy. We can certainly see 10 decimal places of accuracy between what Gauss-Newton gives us with what MATLAB gives us built-in.

Gauss-Newton Solver: Improper assignment with rectangular empty matrix

1 Answers1

Example use: