find roots of determinant

Hello, dear forum members. I have the following problem. I solve the equation f(x)=det(A(x))=0. But because of the large size of the matrix

A, its determinant becomes so large that the matlab considers this to be an infinite number. I temporarily solved this problem by dividing my equation by its value in the initial approximation . However, this dramatically slowed down the calculations. Help me find an alternative way to solve this problem.

[Torsten]

How did you solve the problem until now ? Using "fsolve" ?

[Oleksandr Maksimenko]

Yes, I use fsolve for finding roots

[Oleksandr Maksimenko]

Perhaps this problem can be solved by factoring my matrix, so as not to waste resources on finding the inverse matrix? Or are there other ways to normalize matrices in matlab?

Answer 1

[John D'Errico]

Sadly, too many people misunderstand DET. They think that just because one of the properties of the determinant is thatit is mathematically zero when a matrix is singular, that this is what it should be used to do. Yes, that is brought on because they were probably taught to do exactly that in school, by teachers who were taught the same thing.

DET can be viewed as a tool. It can do some things acceptably well. But like any tool, some are designed for the purpose, and some have a better design than other tools.

For example, suppose I gave you a nail, and a screwdriver, as well as a hammer. Would you try to drive the nail into a block of wood with the screwdriver? You might succeed, eventually, if the wood was balsa wood. But you would be using the wrong tool.

In the floating point arithmetic of a computer, we are faced with the problem that DET does not perform well when testing if a matrix is singular, even though in the world of mathematics, det(A) is zero in theory when A is singular. DET is just the wrong tool to use, and it matters not if you learned the wrong thing in the past. DET is the wrong tool for this purpose.

The biggest problem comes from the scaling behavior of det. For example, is an identiity matrix singular? And, surely, a matrix will all 2s on the diagonal, so just 2*A is not singular? And what if we divide by 2 instead?

A = eye(1200);
det(A)
ans =
     1
     
det(2*A)
ans =
   Inf

det(A/2)
ans =
     0

Surely det should be able to tell you that none of those matrices are singular. In fact, they are all just about the least singular matrices you can imagine. Even worse, I can easily construct a matrix that is clearly numerically singular, yet det will tell you it has in infinite determinant.

So sorry, but in the world of computational mathematics, DET is just a bad seed.

Instead, you should work with what really identifies a matrix as singular. Best would be to work with the singular values of A, but cond will get you close enough, in a way that does not scale poorly like DET. RANK also tells you if the matrix is singular, but RANK returns a discrete quantized value, so it cannot be used in a solver like fminbnd or fzero.

For example, if you want to solve for x, such that A(x) is a singular matrix, then one solution might be to try to drive cond(A(x)) to infinity. Or, you could try to drive 1/cond(A(x)) to zero. And since Afun is always positive, it might be best to just use fminbnd on it. I suppose you might try to miniimize -cond(A(x)). That should work too.

So now perhaps I have convinced you that using DET is a bad idea for this problem. Are there better ways to solve it? Mathematics may come to your aid. Perhaps. For example, suppose you wanted to solve the problem:

A0 = rand(1000);

Now, construct the related matrix A(u), such that A(u) s identical to A0, except that we can perturb the element of A0(1,1) to have any value we wish. Clearly, this probem could be more complex but this is a simple example. So assume we have

A(1,1) = A0(1,1) + u

and A(i,j) = A0(i,j) otherwise.

What value of u willl make this matrix A(u) singular? Is u unique?

det(A0)
ans =
   Inf
cond(A0)
ans =
   5.6904e+05

So see that while A0 is not in itself singular, det returns inf. While cond tells us A0 is not too terribly close to singular.

We can think of the matrix function A(u) as having a simple derivative matrix, of the form:

dA = zeros(size(A0));
dA(1,1) = 1;

Now just use eig.

[V,D] = eig(A0,dA);
D = diag(D);
format long g
D(1)
ans =
          9.39697748820705

We see that D(1) is the perturbation needed to make A0 singular. And V(:,1) gives us the nullspace.

rank(A0 - D(1)*dA)
ans =
   999
norm((A0 - D(1)*dA)*V(:,1))
ans =
      2.21161733072466e-13

So at least for this nice simple problem, all we needed to do was understand how a generalized eigenvalue can solve the problem. We never needed to compute a determinant at all. In fact, since A(u) is linear, no iterative computation was ever needed. However, we also should also be able to write a similar iterative algorithm for a nonlinear problem. Think of it as sort of a Newton's algorithm, where each iteration step would be based on a linear approximation to the perturbation.

[Zack Fifer]

Thank you for this!

Are there any references to cite for this method, or is this something you invented to answer this question?

[John D'Errico]

As I responded via e-mail, it was created on the spot to solve the problem at hand. I might offer Matrix Computations, by Golub & Van Loan, where a discussion is made of generalized eigenvalues.

[Ductho Le]

@John D'Errico what if that function has multiple solutions? fminbnd only returns one solution, furthermore the value of cond is positive thus we can not check whether the function changes sign. Normally, I use islocalmin to find all local minima and check the sign changing at every value to find solution. Do you have any solution for this case?

Answer 2

[James Tursa]

Since the determinant of a matrix is equal to the product of its eigenvalues, maybe you could work with eig(A(x)) instead and try to drive one or more of the eigenvalues to 0. Working directly with the determinant of a large matrix seems like a numericallly instable way to go about this anyway.

[Oleksandr Maksimenko]

Thank you very much for your answer. And how to determine what kind of eigenvalue must be equal zero? Indeed, only one eigenvalue will equal to zero with the root that I want to find.

[John D'Errico]

Please stop adding new answers. You are not answering the question, merely making comments.

"So, I have tried the approach suggested above. Since I did not know which particular eigenvalue must be equal zero, I have solved the equation product of eigenvalues=0 . And it turned out that this approach is even more complicated and requires about 3 times more time to complete the calculations."