How to convert into double?

Hello,

I wrote the following code:

tfail = [5571.760,5573.742,5654.457,6079.693,6081.927,6172.915,6515.064,6517.515,6617.308,7095.558,7098.298,7209.831,7530.929,7533.885,7654.224,7966.300,7969.472,8098.617,8401.671,8405.059,8543.009,8982.166,8985.843,9135.533,9852.908,9857.017,10024.38,10868.774,10873.387,11061.234];
n=length(tfail);

beta_hat = 4.2915822; 
B_hat = 1861.6186657; 
C_hat = 58.9848692;

syms t B beta C
y(t) = (exp(-B/((heaviside(t)-heaviside(t-2000))*(330)+(heaviside(t-2000)-heaviside(t-3000))*(350)+...
    (heaviside(t-3000)-heaviside(t-14000))*(390))))/C;

ogL=0;

for i=1:n 
    tfail(i);
    I(i) = int(y(t),t,0,tfail(i));
    y_new(i)=subs(y,t,tfail(i));
    logL =logL+log((beta*y_new(i)*(I(i))^(beta-1))*exp(-((I(i))^beta)));
end

p = int(y(t),t,0,14000);
u = beta*log(p);
du_dB = diff(u,B);
du_dbeta = diff(u,beta);
du_dC = diff(u,C);
du_dB_sub = subs(du_dB,{beta,B,C},{beta_hat,B_hat,C_hat});
du_dbeta_sub = subs(du_dbeta,{B,C},{B_hat,C_hat});
du_dC_sub = subs(du_dC,{beta,B,C},{beta_hat,B_hat,C_hat});

v=[beta;B;C]; 
H=hessian(logL,v); 
H_negatv=-1*H;

now I would like to calculate the inverse of H_negatv by using:

H_inverse=inv(H_negatv);

But that doesn´t work. So I tried out:

h = 1\H_negatv.

That´s good so far.

But now I would do sth. like that:

w=subs(h,[beta,B,C],[beta_hat,B_hat,C_hat]);
F_direct = w;

In according to calculate:

Var_B_hat_direct = double(F_direct(2,2));

But I can´t do that in MATLAB.

Does somebody have an idea how to solve that problem?

[John D'Errico]

I assume that you meant to type logL = 0, instead of ogL = 0. That may in fact be why you have trouble. I don't know.

But NO NO NO NO! Then you write:

H_inverse=inv(H_negatv); 
But that doesn´t work. So I tried out:

h = 1\H_negatv.

The above does not compute the inverse of a matrix!!!!!!!!!!!!

A = sym(magic(3))
A =
[ 8, 1, 6]
[ 3, 5, 7]
[ 4, 9, 2]

1\A
ans =
[ 8, 1, 6]
[ 3, 5, 7]
[ 4, 9, 2]

So, NO, it is NOT good so far.

[Walter Roberson]

My tests show that inverting H_negatv is not practical on my system (8 Gb) because inverting it runs out of memory. The length of the hessian is about 4 megabytes.

Switching to using \ instead of inv() is not going to help.

As I am pretty sure I said before at some point, you are not looping substituting different values for the symbols in H_negatv, so what you should be doing is substituting in the numeric values and double() the result before doing the inverse, so that you are taking the inverse of a numeric 3 x 3 instead of a large symbolic matrix.

[Max]

Hey John,

but what´s about that:

"Note: MATLAB® computes X^(-1) and inv(X) in the same manner, and both are subject to the same limitations." see, also: http://de.mathworks.com/help/matlab/ref/inv.html

[John D'Errico]

READ WHAT I SAID! Think about what you see in my comment.

As I showed, 1\A does NOT compute a matrix inverse, unless you think that magic(3) is it own inverse.

In fact, when you write the expression 1\A in MATLAB, what you did was to solve the system

1*X = A

The solution to that problem is NOT the inverse of the matrix A. The solution is A.

So giving me lots of documentation that shows that inv(A) is equivalent to x^-1 is irrelevant. None of those sources claim that you will get a matrix inverse from the form 1\A.

PERIOD.

[Walter Roberson]

You need eye(size(A))\A with square nonsingular A to have the equivalent of inv(A)

[Max]

Dear John,

first of all, you´re right concerning the matlab command: 1\A. That doesn´t calculate the inverse matrix.

But instead of doing Ainv = inv(A), I also can do: A_inv = A^(-1):

A = [ 1 2 3; 2 3 4; 1 2 5]
A =

       1     2     3
       2     3     4
       1     2     5

inv(A)

ans =
-3.5000    2.0000    0.5000
    3.0000   -1.0000   -1.0000
   -0.5000         0    0.5000

A_inv
ans =
   -3.5000    2.0000    0.5000
    3.0000   -1.0000   -1.0000
   -0.5000         0    0.5000

My problem is that I can´t do it in MATLAB because of the "out of memory"!!!

[Walter Roberson]

You are not looping substituting different values for the symbols in H_negatv, so what you should be doing is substituting in the numeric values and double() the result before doing the inverse, so that you are taking the inverse of a numeric 3 x 3 instead of a large symbolic matrix.

Answer 1

[Walter Roberson]

Working notes for me.

rational = @(V) sym(V, 'r');
if ismember( exist('hessian'), [2, 3, 5, 6, 8])   %is it an executable function?
  Hess = @(M,V) hessian(M,V)
else
  Hess = @(M,V) maple('Student[VectorCalculus][Hessian]', M, maple('convert', V, 'list'));
end

tfail = rational([5571.760, 5573.742, 5654.457, 6079.693, 6081.927, 6172.915, 6515.064, 6517.515, 6617.308, 7095.558, 7098.298, 7209.831, 7530.929, 7533.885, 7654.224, 7966.300, 7969.472, 8098.617, 8401.671, 8405.059, 8543.009, 8982.166, 8985.843, 9135.533, 9852.908, 9857.017, 10024.38, 10868.774, 10873.387, 11061.234]);

n = length(tfail);

beta_hat = rational(4.2915822); 
B_hat = rational(1861.6186657); 
C_hat = rational(58.9848692);

syms t B beta C
y = (exp(-B/((heaviside(t) - heaviside(t-2000)) * (330) + (heaviside(t-2000) - heaviside(t-3000)) * (350) + (heaviside(t-3000) - heaviside(t-14000)) * (390))))/C;

Z = rational(0);
LogL = Z;

I = rational(zeros(1,n));
y_new = rational(zeros(1,n));
new_term = rational(zeros(1,n));

for i = 1:n 
    I(i) = simplify( int(y, t, Z, tfail(i)) );
    y_new(i) = simplify( subs(y, t, tfail(i)) );
    new_term(i) = log((beta * y_new(i) * (I(i))^(beta-1)) * exp(-((I(i))^beta)));
    logL = logL + new_term(i);
end

p = int(y, t, Z, rational(14000));
u = beta * log(p);
du_dB = diff(u, B);
du_dbeta = diff(u, beta);
du_dC = diff(u, C);
du_dB_sub = subs(du_dB, {beta, B, C}, {beta_hat, B_hat, C_hat});
du_dbeta_sub = subs(du_dbeta, {B,C}, {B_hat,C_hat});
du_dC_sub = subs(du_dC, {beta,B,C}, {beta_hat,B_hat,C_hat});

v = [beta; B; C]; 
H = Hess(logL, v);

%the next will probably fail, running out of memory
Hs = simplify(H);
H_negatv = -1*Hs;

[Walter Roberson]

The above is not intended as a solution: it is a recoding for compatibility with the Maple interface. On my system it will run out of memory attempting to simplify H. And if you skip the simplify() step and go on to taking the inv() it will run out of memory trying to take the inverse.

Your expression is simply too large to take a reasonable symbolic inverse of. But since you have specific numeric values to substitute in, you should substitute those specific numeric values into H, take double() of the result (since it will be completely numeric) and inv() that purely numeric interface.

Except, of course, we do not recommend that you take inv() for any reason other than proving that you can.