MATLAB: How to solve the determinant of large-sized symbolic matrix: det(A)=0

determinantsolvesymbolic matrix

I have two n*n-sized matrix, Y is a constant tri-diagonal matrix,

where all the elements in Y is known and constant.

D is a diagoanl matrix,

where only d is unknown, the other numbers are known (e_n, f_n are known)

I want to find the complex value of d that can satisfy det (D-Y)=0.

Here is what I did:

I define a matrix Y (20*20)
I define d is a symbolic value: syms d
I define D is a symbolic matrix : sym('D', [20 20]), with symbolic value d
I define K=det(D-Y);
I solve p=double(solve(K,d))

It can be solved if the matrix size is small, for example 5*5, but if the matrix size is large, for example 20*20, it is super slow.

I guess the reason is that, there is square roots expressions in matrix D which makes the determinant K cannot be expressed as a polymonid of d.

Is there any good method to solve fastly? Thanks for help!

Best Answer

First, you don't want to solve for a zero determinant! Period. This will never be a well-posed problem, simply because we can trivially come up with a variety of matrices that have an essentially zero determinant yet are not singular, or others that have an essentially infinite determinant, yet they are numerically singular. For example, consider this simple 20x20 diagonal matrix.

det(diag([repmat(1e10,1,19),1]))
ans = 1.0000e+190

It is not singular, yet it has a determinant of 1e190. This next, similarly constructed matrix is also not singular.

det(diag([repmat(1e-10,1,19),1]))
ans = 1.0000e-190

Yet the determinant is 1e-190. With only slightly more effort, I could have made the determinants underflow to zero in double precision, or overflow to inf.

Yes, I know your teacher told you to use the determinant in some long half-forgotten class. What they either forgot to tell you, (or you just forgot) is not to use it on any real world problem. Determinants are great things for students to play with. But they will always lead you into trouble.

Next, you cannot solve this symbolically. Sorry, but you cannot. The answer most of the time when a new user does not know the value of some variable, is they make it into a symbolic problem. It looks impressive, but that does not mean this is a viable way to solve the problem.

Once you make it too high in dimension, the symbolic expressions get immensely large. Instead, you need to recognize this is a simple 2 dimensional problem, that can be solved simply using double precision arithmetic. It is a 2-d problem because the complex variable d is really 2 variables, comprising the real and an imaginary parts of d. Anyway, since you cannot find a symbolic solution, meaning solve is useless here for you, you need to recognize that no algebraic solution exists. You could never do that, because you really can formulate this as a super high degree polynomial, and once you get past degree 4, Abel-Ruffini proved the impossibility of this long ago.

So, how would I approach this problem? I lack your numbers. So I will make some up.

Y = tril(triu(randn(20),-1),1);
E = randn(20,1);
F = randn(20,1);

So Y is some randomly garbage tridiabonal matrix. E and F are randomly garbage vectors, each of length 20. Should some of them be complex? Yeah, maybe. I don't know, and you don't seem to tell me. Honestly, I don't care, because that is not relevant to my solution.

Now, I need to create the martrix D, as a function of d. It is not vectorized, because it returns an entire matrix. (I could have done it more elegantly I suppose, but that would have taken more work to avoid a warning message.) As a function of two unknowns, the real and imaginary parts of d,

Dmat = @(dr,di) diag(sqrt((dr + i*di).^2 - E)./F);

Next, is the resulting matrix D-Y singular? As an example, I'll pick some values of the two unknowns. The matrix will be considered to be numerically singular in double precision, if rcond((D-Y) is less than eps, or at least, on the order of eps.

rcond(Dmat(1,3) - Y)
ans = 0.0077

In fact, we can plot the reciprocal condition number of D-Y, as a function of the two variables.

fsurf(@(dr,di) rcond(Dmat(dr,di) - Y))
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.
xlabel real(d)
ylabel imag(d)
zlabel 'Reciprocal condition number'
●

So just use a minimizer to find those locations where rcond approaches zero. Even a high quality optimizer (sarcasm inserted) like fminsearch would suffice. No matter what you do, you will need to search for each location that produces a singularity, one at a time. Here is one solution on this problem:

[dri,rcondval,exitflag] = fminsearch(@(dri) rcond(Dmat(dri(1),dri(2)) - Y),[1 1])
dri = 1×2
    1.4278    0.8799
●
rcondval = 2.1118e-07
exitflag = 1

Choose different starting points for the optimizer and you will find other solutions. It looks like there are many.

Again, better code written on my part would even figure out how to vectorize it to avoid the warning messages, but it is just not significant here.

Essentially, this is a simple problem to solve, IF you approach it the correct way. Or you can use det and syms, and get crap that you will never be able to use/trust, even if it does eventually return a result next year. Your choice.