I don't understand how you want to define the hessian with respect to "par" for vectors x and x0.
Do you want to create 100 x 100 Hessians of size 2x2, one for each pair (x(i),x0(j)) ?
An issue with matlabFunction
5 views (last 30 days)
Show older comments
Hello,
I have a question and I think the best is to explain via an example. So, consider the following
mu=@(x,par)par(1).*x;sigma=@(x,par)par(2);
syms x0 x dt
par = sym('par', [1, 2]);
LL = -1 / 2 * (log(2 * pi * dt * sigma(x0, par)^2) + ((x - x0 - mu(x0, par) * dt) / (sigma(x0, par) * sqrt(dt)))^2);
h_LL = hessian(LL,par);
h_LL=matlabFunction(h_LL)
Now, assume that V is a 1-by-n given vector and x0=V(1:n-1), x=V(2:n), dt and par's are given scalars. for instance, assume that dt=1, x0=1:100, x=2:101, par1=par2=1. What I want to do is to calculate the following sum over the elements of V:
h_LL(1,1,1,1,2)+h_LL(1,1,1,2,3)+h_LL(1,1,1,3,4)+ ... + h_LL(1,1,1,100,101)
Of course I can calculate this sum using a for-loop. But, my actual 'n' and the number of my actual parameters are very big making this scheme impractical. The idea that comes to my mind is to find a way to remove the 'reshape' in the above function handle h_LL and supply the entire vevtors x and x0 rather than doing so element by element.
Unfortunately, I haver no control over matlabFunction and cannot tell it please give me a matrix rather than applying reshape to a scalar.
Any idea?
Thanks in advance,
Babak
2 Comments
Accepted Answer
Dyuman Joshi
on 11 Dec 2023
Edited: Dyuman Joshi
on 11 Dec 2023
mu=@(x,par)par(1).*x;
sigma=@(x,par)par(2);
syms x0 x dt
par = sym('par', [1, 2]);
LL = -1 / 2 * (log(2 * pi * dt * sigma(x0, par)^2) + ((x - x0 - mu(x0, par) * dt) / (sigma(x0, par) * sqrt(dt)))^2);
temp = hessian(LL, par)
%Reshape the hessian as a column vector
h_LL = reshape(temp, [], 1)
%Convert to a function handle
fun=matlabFunction(h_LL);
%Sample inputs
DT=1; X0=1:100; X=2:101; PAR1=1; PAR2=1;
%Provide x0 and x as row vectors
out = fun(DT, PAR1, PAR2, X, X0)
%Convert the output to a 3D array,
%where each 2D matrix is a corresponding output
out = reshape(out, 2, 2, [])
%For comparison
subs(temp, [dt x0 x par], [DT X0(1) X(1) PAR1 PAR2])
subs(temp, [dt x0 x par], [DT X0(2) X(2) PAR1 PAR2])
2 Comments
More Answers (1)
Walter Roberson
on 11 Dec 2023
matlabFunction can produce vectorized code only if all of the following are true:
- the expression is scalar for scalar inputs
- the expression does not use piecewise()
- the expression does not involve an integral, unless the integral is one that has a closed-form solution. This applies regardless of whether the symbolic parameters appear in the expression to be integrated, or in the bounds of the integral (but for different reasons for those two situations)
- the expression does not use symmatrix (matlabFunction refuses to process those)
I might be forgetting some other cases.
The idea that comes to my mind is to find a way to remove the 'reshape' in the above function handle h_LL
You would have to do that by post processing, such as by having matlabFunction use the 'File' option and edit the file afterwards. Unless you wanted to risk putting a reshape.m on your MATLAB path, but that would be distinctly risky.
So what is the official solution for this kind of problem?
- You can loop; or
- You can use arrayfun() to apply the function to "corresponding" matrix inputs, which will automatically wrap the results; you would 'uniform', 0 to get a cell array of results. You could reshape() or permute() that cell to be 1 x 1 x n and then use cell2mat() -- or you could do cat(3, ThatCell{:}) to concatenate along the third dimension.
0 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)