function [g coeff]=heat_hermite(f,L,sigma,m)
%--------------------------------------------------------------------------
% Heat diffusion using Hermite polynomial expansion
%
% f :       input data of size N x d (N: data size, d: number of data)
% L :       Laplace-Beltrami operator (N x N sparse matrix)
% sigma :   diffusion time 
% m :       expansion degree
% g :       heat diffusion of each column of f
% coeff :   expansion coefficients
%
%
% The mathematical detail is given in [1].
% The recurrence relation part is modified from the Spectral Graph Wavelet 
% Transform (SGWT) [2] toolbox (https://wiki.epfl.ch/sgwt).
%
%
% References:
% [1] Huang, S.-G., Lyu, I., Qiu, A., Chung, M.K.: 
% Fast polynomial approximation to heat diffusion in manifolds, 
% 2019 Medical Image Computing and Computer Assisted Intervention (MICCAI)
% http://pages.stat.wisc.edu/~mchung/papers/huang.2019.MICCAI.pdf
%
% [2] Hammond, D.K., Vandergheynst, P., Gribonval, R.: Wavelets on graphs via 
% spectral graph theory. Applied and Computational Harmonic Analysis 30, 129-150, 2011
%
% 
% (C) 2019 Shih-Gu Huang shihgu@gmail.com
%          Moo K. Chung  mkchung@wisc.edu
%          Universtiy of Wisconsin-Madison
%
% Update history:
%     July 29, 2019 created by Huang
%--------------------------------------------------------------------------


% expansion coefficients c_{sigma,n} for n=0,1,...,m
n=[0:m].';
coeff=zeros(m+1,1);
coeff(1)=exp(sigma^2/4);
for n=1:m
coeff(n+1)=coeff(n)*(-sigma/2)/n;
end

% initial conditions of recurrence relation 
Pf_old =0;    % P_{-1}(L)f = 0
Pf =f;       % P_{0}(L)f  = f
g=0;


for n=0:m
    % heat diffusion approximation
    g =g +coeff(n+1)*Pf;    
    % recurrence relation: P_{n+1}(L)f =2*L*P_{n}(L)f -2n P_{n-1}(L)f    
    Pf_new =2*L*Pf -2*n*Pf_old;   
    Pf_old =Pf;
    Pf =Pf_new;
end

end