Unstable Matrix exponential when solving ODE

d_values = 10:110;  % Range of d values to scan
t = 0.6;          % Time value
first_element = zeros(length(d_values),1);  % Preallocate array to store the first element of vec1for idx = 1:length(d_values)
for idx =1:length(d_values)    
d = d_values(idx);
    A = diag(sqrt(1:d-1), 1);
    Ad = A';
    A3 = A^3;
    Ad3 = A3';
    H = (Ad3 + A3); 
    vec0 = zeros(d,1);
    vec0(1) = 1;
    U_tsq = expm(-1i * t * H);
    vec1 = U_tsq * vec0;
    first_element(idx) = vec1(1);
end
% Plot d vs the first element of vec1
plot(d_values, abs(first_element), '-o');
xlabel('matrix size');
ylim([0,1]);
ylabel('Abs of the first element of vec1');

grid on;

The above should calculate the matrix exponential which is the propagator of the Hamiltonian where are the annihilation and creation operators of the harmonic oscillator. For various system size truncations d, I get different results for the U|0> transition element where |0> = (1,0,...,0) (i.e. numerical instability). Suggestions on tackling this problem would be helpful.

Answer 1

[Ayush]

Hi @ibes,

Numerical instability can be a common issue when calculating the matrix exponential, especially for large matrices or when dealing with complex numbers.

In your code, the high condition number suggests that the matrix is ill-conditioned, which can lead to such instability. The condition number can be evaluated using MATLAB's “cond” function. You can read more about “cond” function here: https://www.mathworks.com/help/matlab/ref/cond.html

To mitigate numerical instability, you can consider the following strategies:

1. Matrix Balancing: You can use matrix balancing to improve the conditioning of the matrix before computing the exponential. MATLAB's “balance” function can help with this. You can read more about it here: https://www.mathworks.com/help/matlab/ref/balance.html
2. Regularization: Consider adding a small regularization term to the matrix to improve its condition number. This can be done by adding a small multiple of the identity matrix, “epsilon * eye(d)”, to H.
3. Alternative Matrix Exponential Methods: You can also use “Krylov subspace” techniques or other specialized libraries designed for stable matrix exponential calculations. You can read about it here: https://www.mathworks.com/matlabcentral/fileexchange/45170-jacobian-free-newton-krylov-jfnk-method

Here’s my modified version of your code involving regularization and matrix balancing to reduce condition number and bring numerical stability.

d_values = 10:110;  % Range of d values to scan 
t = 0.6;            % Time value 
first_element = zeros(length(d_values),1);  % Preallocate array to store the first element of vec1 

for idx = 1:length(d_values) 
    d = d_values(idx); 
    A = diag(sqrt(1:d-1), 1); 
    Ad = A'; 
    A3 = A^3; 
    Ad3 = A3'; 
    H = (Ad3 + A3);  
    
    % Regularization (small identity matrix addition) 
    epsilon = 1e-10;  % Regularization parameter, adjust as needed 
    H_reg = H + epsilon * eye(d); 
    
    % Balance the matrix 
    [H_bal, D] = balance(H_reg); 
    
    % Check condition number 
    cond_number = cond(H_bal); 
    fprintf('Condition number for d = %d: %e\n', d, cond_number); 
    
    vec0 = zeros(d,1); 
    vec0(1) = 1; 
    
    % Compute the matrix exponential 
    try 
        U_tsq = expm(-1i * t * H_bal); 
    catch ME 
        warning('Matrix exponential computation failed for d = %d: %s', d, ME.message); 
        continue; 
    end 
    
    vec1 = U_tsq * vec0; 
    first_element(idx) = vec1(1); 
end 

% Plot d vs the first element of vec1 
plot(d_values, abs(first_element), '-o'); 
xlabel('matrix size'); 
ylim([0,1]); 
ylabel('Abs of the first element of vec1'); 
grid on; 

Hope it helps!

[ibes]

Thank you for your suggestion. I tried your solution. But as far as i understand it doesnt give the correct result. I tried to fix it: The balancing arguments are flipped. It should be

[D, H_bal] = balance(H_reg)

Further

U_tsq = D exp(-i t H_bal)/D

But then the numerical instabilities are back.