0
$\begingroup$

I have an expression (let's say determinant of matrix A) expressed in symbolic form in terms of 2 decision variables x1, x2 and 2 parameters q1 and q2. I'm minimizing this using fmincon for different start point and parameter value combinations.

For certain combinations of x1, x2, q1, q2 the matrix A is not invertible so I would want determinant evaluation to be zero but Matlab just gives a very small number like 9.7166e-91.

If I had created the matrix A with numeric values for which inverse doesnt exist and then take the inverse it does give me the usual warning "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate." This warning is very good as I know that the output can be discarded.

But I need to work with symbolic expression as this is what needs to be fed to fmincon.

So what approach should I use in order to know that certain outputs from the numeric evaluation of a symbolic expression should be discarded?Like am happy to have "NaN" when inverse doesn't exist. Hope there is a trick which doesnt slow down my program

Note, the Digits setting is 32 so should I use this and infer that if the determinant result is between -1E-32 and 1E-32 then result should be discarded ? (this approach does not seem correct)

Any suggestions are welcome.

%Providing code below in order to replicate easily;

syms x1 x2 q1 q2;% x1 and x2 are the decision variables while q1 and q2 are the parameters. For different values of parameters I'm trying to minimize a function using fmincon


%matrix A is DesignMatrixMult. Example provided below for illustrative purposes but its form can change depending on user input;
 DesignMatrixMult = [                    ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 - (exp(-q1*x1) - exp(-q2*x1))/(q1 - q2) + (q1*x1*exp(-q1*x1))/(q1 - q2))^2 + ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 - (exp(-q1*x2) - exp(-q2*x2))/(q1 - q2) + (q1*x2*exp(-q1*x2))/(q1 - q2))^2, - ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 + (q1*x1*exp(-q2*x1))/(q1 - q2))*((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 - (exp(-q1*x1) - exp(-q2*x1))/(q1 - q2) + (q1*x1*exp(-q1*x1))/(q1 - q2)) - ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1
    - q2)^2 + (q1*x2*exp(-q2*x2))/(q1 - q2))*((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 - (exp(-q1*x2) - exp(-q2*x2))/(q1 - q2) + (q1*x2*exp(-q1*x2))/(q1 - q2));
     - ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 + (q1*x1*exp(-q2*x1))/(q1 - q2))*((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 - (exp(-q1*x1) - exp(-q2*x1))/(q1 - q2) + (q1*x1*exp(-q1*x1))/(q1 - q2)) - ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1
    - q2)^2 + (q1*x2*exp(-q2*x2))/(q1 - q2))*((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 - (exp(-q1*x2) - exp(-q2*x2))/(q1 - q2) + (q1*x2*exp(-q1*x2))/(q1 - q2)),                                        ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 + (q1*x1*exp(-q2*x1))/(q1 - q2))^2 + ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 + (q1*x2*exp(-q2*x2))/(q1 - q2))^2];


% determinant is calculated below;  
DCritn = simplify(simplify(det(DesignMatrixMult),'IgnoreAnalyticConstraints', true,'steps',500),'full');
%creating an anonymous function handle to DCritn;
DCritnF = matlabFunction(DCritn,'vars',{x1,x2,q1,q2},'file','');



% few lines of code for generating different possible values of the parameter q1 and q2;
% As an example, one set of values is provided below;
q1=0.9287; q2 = 0.83;

For the above parameter values, the function DCritnF is minimized using fmincon;

The challenge I have is that for certain combinations of x1 and x2, DCritnF gives non-zero result even though I know that the inverse of DesignMatrixMult (matrix A) does not exist. For example, DCritnF is called with value of x1 = 42.5 and x2 = 95, it returns 9.7166e-91

On the other hand, if I do the following,

TrialCheck = double(subs(DesignMatrixMult));
TrialCheckInverse = inv(TrialCheck);

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.404300e-18. ans =

1.0e+42 *

5.4788 -1.3938

-1.3938 0.3546

Also one smaller question;if I had done

TrialCheck_Alternate = double(subs(DesignMatrixMult))
TrialCheckDeterm_Alternate = det(TrialCheck_Alternate)

Then I get no warning? If for inverse of matrix it issues a warning then shouldn't it also issue warning while finding determinant of that matrix?

TrialCheck_Alternate =

1.0e-25 *

0.0139 0.0546 0.0546 0.2146

TrialCheckDeterm_Alternate =

3.9175e-69

$\endgroup$
  • $\begingroup$ I was searching the web and it seems one way could be to use eps and if true I would need to compare the absolute value of determinant with eps? $\endgroup$ – Hari Jul 11 '14 at 18:25
0
$\begingroup$

I'm guessing what you want is some way to tell that your square matrix $A$ is either singular or non-singular. For that purpose, I'd use the rank function in MATLAB, which is internally using the SVD with some thresholding to determine the number of nonzero singular values.

$\endgroup$
  • $\begingroup$ Thanks for the response. The challenge is that I need to optimize (actually maximize and not minimize) the determinant of A (which is DCritnF) and so I feed negative of DCritnF to fmincon. If I have to check rank then I cannot just use DCritnF to evaluate the function but also construct the matrix A using the present values of x and q and then apply rank function (increased execution time). $\endgroup$ – Hari Jul 11 '14 at 20:34
  • $\begingroup$ The problem with just comparing the value of the determinant with machine precision is that it's not incredibly meaningful. For instance, a scaled identity matrix could have a small determinant, but should yield a numerically well-behaved inverse. It's not clear to me how you get what you want without your matrix unless you're willing to accept potentially poor results. $\endgroup$ – Geoff Oxberry Jul 12 '14 at 0:14
  • $\begingroup$ I understand now that rank is the right way to approach this. Thanks for the clarification. I will try to work on the right way to integrate checking of rank prior to evaluation of objective function. Will get back in case of further questions. $\endgroup$ – Hari Jul 12 '14 at 1:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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