Instantly share code, notes, and snippets.

# higham/funm_randomized

Last active June 8, 2018 13:46
Show Gist options
• Save higham/6c00f62e48c1b0116f2e9a8f43f2e02a to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 function F = funm_randomized(A,fun) %FUNM_RANDOMIZED Evaluate general matrix function using % randomized approximate diagonalization method of Davies (2007). tol = 8*eps(A(1,1)); % Tolerance for single or double precision A. E = randn(size(A)); [V,D] = eig(A + (tol*norm(A,'fro')/norm(E,'fro'))*E); F = V*diag(fun(diag(D)))/V;

### advanpix commented Apr 7, 2017

The code is easy to convert to work with any precision (not only single or double).
The only one change is to create random matrix of the same type as A:

`E = randn(size(A),class(A));`

Now our code has became precision-independent and can be used with double precision:

``````>> A = rand(5,'double');
>> norm(funm_randomized(A,@exp)-expm(A),1)/norm(A,1)
ans =
7.75329562589415e-15
``````

``````>> mp.Digits(34);
>> A = rand(5,'mp');
>> norm(funm_randomized(A,@exp)-expm(A),1)/norm(A,1)
ans =
7.578684267133064009245947356284363e-33
``````

or any other:

``````>> mp.Digits(1000);
>> A = rand(5,'mp');
>> norm(funm_randomized(A,@exp)-expm(A),1)/norm(A,1)
ans =
3.16466380.....10447022701e-998
``````

### higham commented Apr 7, 2017

Thanks, Pavel - that's an excellent point. E should indeed be of the same precision as A and your change achieves that very neatly.