MATLAB: Basic iterative schemes in matlab

linear

Hello,

I am becoming acquainted with solving basic linear problems in matlab iteratively. So, I apologize for this very basic question.

I am writing a program to solve Ax=f, for a randomly generated matrix (say, of size 6).

I want to extract to break A into M and N, where A=M-N.

I have written a code for this, and I would like to determine the number of iterations that are required before the solution converges. However, my solution does not converge! It blows up to infinity. I have checked to make sure the matrix is non-singular, and I think I have the iterative scheme correct. I would be very grateful if someone could please check my work and let me know where I might be going wrong. Here is my code:

scale=6; 
%A is a random matrix 
A = randn(scale,scale);
if det(A)==0
    fprintf('Matrix is singular');
    %break
else
    fprintf('Matrix is non singular');
end
n=size(A,1);
%need to construct M and subtract N from it
%M contains diagonal elements, here a tridiagonal;
%N contains off diagonal elements of A
M=tril(triu(A,-1),1); %tri diagonal
%tril(triu(A,-2),2) %penta diagonal
N=-(M-A);
%B is a random name for a check to ensure that A=M-N
B=M-N; %Check that we infact get back matrix A
%Solve A*x=f
f=randn(length(M),1);
x=zeros(length(M),1); %starting vector
k=1; %iteration parameter
r=f-A*x(:,k); 
tol=1e-3; %tolernace of the convergence, dictated by the 2-norm of r
while norm(r)>=tol
    x(:,k+1)=M\-N*(x(:,k) + f)
    r(:,k)=f-A*x(:,k);
    norm(r)
    k=k+1;
end

Edit:

Added a condition:

if abs(max(E))>=1
    fprintf('Conditions do not hold')
    break
end

However, there are times when the random matrix A satisfied this condition, and runs through the while look, and still blows up! Any ideas?

Best Answer

(Note: I've not checked to see if your iterative code actually represents the scheme you seem to want to use. Really, I had no need to do so, because the answer is, it WILL diverge almost always for such a random mnatrix.)

Your question comes down to, under what circumstances would such a scheme diverge? Perhaps you want to consider if the Jacobi method always converges, for any matrix. (No, it will not in general converge for such a random matrix.)

As the wiki page points out,

scale=6; 
A = randn(scale,scale);
M=tril(triu(A,-1),1); %tri diagonal
N=-(M-A);

Now, what does that page point out as a general requirement for convergence for such a scheme?

max(abs(eig(M\N)))
ans =
   3.3047

Significantly greater than 1. Not gonna converge.

Related Solutions

MATLAB: Is the rank function reliable

Many things to discuss here.

1. Yes, rank is quite reliable. Is it perfect? Well, no numerical method can be perfect. One can trivially generate a matrix that is not singular, but in double precision determining the rank using any simple scheme will fail to determine the singular nature properly.

A = diag([1 1e-20])
A =
            1            0
            0        1e-20
rank(A)
ans =
     1

So A is clearly diagonal, with non-zero elements, but in double precision, it is numerically singular as far as rank can tell.

However, rank IS essentially as good a scheme as you can get in double precision. If rank tells you the matrix is singular in double precision, then effectively it is so. Anyway, floating point arithmetic limits what you can know about any such matrix.

2. Can you use cond instead of rank? Of course you can, but if rank tells you the matrix has full rank, then cond would do so too. The extra thing that cond tells you is how CLOSE the matrix is to singular, but you still need to make a judgment about the singularity. Similarly, looking at the singular values of the matrix would tell you something, IF you understand what you are looking at.

Regardless, just use rank. For the other tools that require judgement about how close a number is to big, if you don't know what that number should be, then you are FAR better off to just use rank! (My brother-in-law once asked me what a certain tool did, and if he needed it in his shop. My response was essentially that if he had no idea what the tool was, he did not need it.)

3. Evaluating det(A'*A), and wondering how close it need to be to zero is INSANE. The determinant is NEVER a good test of singularity. This is a myth perpetrated on students by authors of textbooks and papers, authors who sadly don't know any better. They keep reading papers by other authors of the same ilk, and repeating what they see. Of course, this is an infinite loop.

The determinant is NEVER a good test of singularity. PERIOD. PERIOD squared.

Once can trivially create a matrix that has any determinant you wish, simply by scaling that matrix by any arbitrary constant. Clearly scaling a matrix by a constant has NO effect on the rank of that matrix, however it DOES strongly affect the determinant. In fact, for an NxN matrix A

det(A) == 1

then

det(k*A) == k^N

So if k is small, then det(k*A) will be tiny! Likewise, if k is large, then det(k*A) will be immense. Since the value of k here has no impact on the singularity of k*A, clearly the det is meaningless to know if the matrix is singular. PERIOD. I don't care how many papers you have read by authors who have no understanding of the issues, how many books. There are a lot of fools (sorry, but true) out there who will endlessly repeat what they have read, not understanding the mathematics.

4. Of course, computing det(A'*A) is worse yet. Here we are squaring the condition of the matrix, so if we could not determine if A was singular using det, why in the name of god and little green apples would forming A'*A help?

5. Similarly, comparing rank(A) to rank(A'*A) is meaningless.

A = diag([1 1e-10]);
rank(A)
ans =
     2
rank(A'*A)
ans =
     1

Was A singular? Of course not. Rank told you so, but you can screw things up by forming A'*A when there was no need to do so. In fact, it is almost NEVER a good idea to form A'*A. Almost any problem that "needs" A'*A will almost always be better served if you read some of your old linear algebra 101 notes, in terms of a factorization of A.

DON'T compute A'*A. Instead, you might compute the QR factorization of A, then use that to form what is essentially a Cholesky factorization of A'*A, BUT keep it in factored form! I.e., if

A = Q*R

then

A'*A = R'*Q'*Q*R = R'*R

Again, DON'T compute R'*R. That would just get you in the same hot water as before. Keep R'*R in factored form. Anything you need to do with that can be gotten from R, and more efficiently AND more accurately. (READ YOUR NOTES! If you never took linear algebra, it is high time to do so, or at least read a book on the subject if you are asking questions like this.

Sorry about the rant. It is not directed at you personally. But if nobody ever tells you the issues and the reasons to avoid things like det, etc., you will never learn, and you will likely perpetuate the myth. There are GOOD numerical schemes to be found for many of the things we do. At the same time, there are many schemes out there that are numerically useless or worse, because they are dangerous if people use them without understanding the issues.

set('soapbox','off','rant','off')

MATLAB: Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

What do you mean by "tolerance"?

A QR solution has no need for a tolerance, nor does linsolve, as they are not iterative methods.

Anyway, you cannot specify the desired accuracy of a solution in a linear system like this, solved using QR. Just wanting a more accurate solution is not relevant. The solution exists, and is unique, or there are infinitely many solutions, all of which are equally good (or poor), or there is no exact solution, and that which qr/linsolve returns is entirely valid.

Perhaps you have a singular system, or one that is nearly singular. In that case, a tolerance is still not really relevant. Yes, you could decide which singular values to discard, in case of a pinv style solution. But that is not a "tolerance" in the way tolerances are traditionally used. And if your system is singular, then again, tolerances have no real meaning.

As far as your matrix being non-square, that just means (assuming you have an over-determined system) that you will never be able to achieve an exact solution. (Well, not entirely true. It IS possible that an over-determined system has an exact solution. But that is terribly rare.) So if no exact solution exists, then specifying a tolerance is again a meaningless goal.

If your problem is non-square in the sense that it is under-determined, then again, a tolerance is irrelevant. A solution may exist, or not. If one solution does exist, then there will be infinitely many solutions, all equally good. A tolerance is meaningless here.

I think that perhaps you are being confused, thinking that the solution to a linear system is like that of a nonlinear system, where you will employ convergence tolerances, etc. Or, perhaps you think that simply by cranking down on a tolerance, you can always gain an arbitrarily good solution. (Many people seem to think that. They are flat out wrong.) Or perhaps you are thinking that a solution like that provided by linsolve/qr is like an iterative solver, like that from lsqr.

If this has not cleared things up, then I will suggest you need to explain CLEARLY what you want to do, and why you think that a tolerance is a meaningful thing to ask for here. As well, you need to explain why it is that you think what I have said above does not apply.

Best Answer

Related Solutions

MATLAB: Is the rank function reliable

MATLAB: Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

Related Question