Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I have made following matlab code for computing the Moore-Penrose inverse of a given matrix A.

A = randn(100)            % given matrix

beta = 1/norm(A,2)^2

x0 = beta.*A'             % initial approximation, A' is the transpose of matrix A

k = 0;

iter = 0

f = 1 ;

I = eye(100);


while (f > 1.0e-007)

x1 = 2*x0 - x0*A*x0;   % x1 is approximation of Moore-Penrose inverse of matrix A

iter = iter+1

a = norm(A*x1*A-A, 2);  %  error norm

 b = norm(x1*A*x1 - x1);`   % error norm

c = norm ((A*x1)' - A*x1, 2 );   % error norm

d = norm ((x1*A)' - x1 *A, 2);   % error norm

B = [a, b, c, d];

f = max(B);

x0 = x1;

end

Since matrix A is changing after every iterations hence values of a,b ,c and d are also changing. I am trying to compute average of values a, b, c and d after hundred repetitions but I am unable to do so. Could anybody help me? I would be very much thankful to you.

share|improve this question
    
Change your code to "B(iter,:) = [a,b,c,d];" and "f = max(B(iter,:));" This will store the value of a,b,c,d for each iteration as successive rows of the matrix B. After the loop ends, you can find the averages using "mean_a = mean(B(:,1));" and similarly for the others. –  in_wolfram_we_trust Sep 25 '12 at 10:37
    
@in_wolfram_we_trust i did this although value of a is in the power s of $10^{-13}$ but mean a is giving me results like 1.33...while it should be similar to a..since it is simply mean of values of a after 100 repetitions... –  srijan Sep 25 '12 at 10:50
1  
I've just run the code as well at found mean_a = 0.8126. The difference in values is due to a different initial seed, but both these answers are acceptable. If you use plot(B(:,1)) it shows you how the value of a has changed - while the value of a at the last iteration is very small (in my case ~10^-11) it starts off quite large (in my case ~10). Averaging over the whole set of values gives the answer you think is too big. –  in_wolfram_we_trust Sep 25 '12 at 11:04
1  
@srijan wait you mean after 100 iterations or after rerunning the code with 100 different matrices A ? –  sonystarmap Apr 11 '13 at 11:32
1  
So lets say for a single matrix your method converges, you get $a=10^{-8}$ and for the next matrix your method converges with $a=10^{-9}$ you would want to get $mean(a) = 5 \cdot 10^{-8}$ ? –  sonystarmap Apr 11 '13 at 11:45
show 8 more comments

1 Answer

up vote 1 down vote accepted
+50

I think the reason, why the mentioned approaches didn't work was because everyone was assuming ( me too at first ) that you want to have the mean value for a during the iterations with one single matrix. However, you seem to want the mean value of a at the end of the algorithm, averaged for several runs with different matrices.

If I understand this correct, the below code should work. (Note $M=10$ instead of $M=100$, line 2)

clc              % cleans the output screen, I don't like old information...
M = 10;          % number of matrices for which perform the calculation
E = zeros(M,4);  % we will store for all M iterations all the 4 values

fprintf('\nPerforming %i iterations in total.\n',M);  % this is just output..

%now the important part, we will perform your calculation M times 
%therefore, we loop over your calculation with a for loop
for i = 1:M

    fprintf('Iteration %i is running...\n' ,i ) % output..

    A = rand(100);            % given matrix
    beta = 1/norm(A,2)^2;
    x0 = beta.*A';             % initial approximation, A' is the transpose of matrix A
    k = 0;
    iter = 0;
    f = 1 ;
    I = eye(100);


    while (f > 1.0e-7)

        x1 = 2*x0 - x0*A*x0;   % x1 is approximation of Moore-Penrose inverse of matrix A
        iter = iter+1;
        a = norm(A*x1*A-A, 2);  %  error norm
        b = norm(x1*A*x1 - x1);   % error norm
        c = norm ((A*x1)' - A*x1, 2 );   % error norm
        d = norm ((x1*A)' - x1 *A, 2);   % error norm

        B = [a, b, c, d];


        f = max(B);

        x0 = x1;

    end

    % at this point, your calculation has been performed for a matrix A.
    % We will now store  the values for a,b,c,d of the last iteration in the matrix E
    E(i,:) = B; 
end

% here we print the mean value of a, which is mean(E(:,1))
fprintf('\nMean value of a =%s\n', mean(E(:,1)));
share|improve this answer
    
I have no words to thank you. I just want your little explanation about the code. I am not very much good in matlab coding. Thanks again. :) –  srijan Apr 11 '13 at 12:05
1  
No hurry. I will add some comments to the code soon. –  sonystarmap Apr 11 '13 at 12:14
1  
Thanks again sir. –  srijan Apr 11 '13 at 12:19
1  
@srijan Comments added. I think that everybody else also could have helped, but the problem was that we first missinterpreted your question ;) –  sonystarmap Apr 11 '13 at 12:20
1  
Dear sir everything is clear now. No queries. Heartily thanks. :) –  srijan Apr 11 '13 at 12:28
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.