I am trying to code up a numerical algorithm which takes a spectral function of the form $$c(\zeta) = w_0 +\sum_{m=1}^N \frac{w_m}{\lambda_m+\zeta}$$ into a continued fraction of the form $$c(\zeta) = d_0+\frac{1}{a_0\zeta+}\;\;\frac{1}{b_1+}\;\;\frac{1}{a_1\zeta+}\;\;\frac{1}{b_2+}\;\;\frac{1}{...}$$ I am using a paper by Weidelt (2005) where he details the numerical algorithm he uses (in Section 2.2), but I can still not make heads or tails of his algorithm. From the spectral function, he builds a rational polynomial with coefficients $p_n$ and from there he is able to translate the $p_n$ coefficients into $a_n$ and $b_n$ coefficients of the continued fraction. The algorithm is as follows:

Let $\beta_0 = 0$, $p_{-1}(\lambda)=0$ and $p_0(\lambda)=1/\sqrt{s_0}$ where $s_0 = \sum_{m=1}^N w_m$.

For $n= 0,1,2,...$ the polynomial can be solved recursively:

i) $\alpha_n = (\lambda p_n, p_n)$

ii) $p^*_{n+1}(\lambda)=(\lambda-\alpha_n)p_n(\lambda)-\beta_np_{n-1}$

iii) $\beta_{n+1}=\sqrt{(p^*_{n+1},p^*_{n+1})}$

iv) $p_{n+1}(\lambda) = p^*_{n+1}/\beta_{n+1}$

After finding the polynomials, $a_n$ and $b_n$ can be solved for easily.

I primarily don't understand 3 things:

1) What do the parentheses in Step i and iii refer to?

2) What is the $n$ subscript referring to?

3) For discrete $\lambda$, how does a vector of $\lambda$ values work into the algorithm? His algorithm seems to suggest $\lambda$ is either a singular number or a continuous variable.

Any help is appreciated and if anyone has the will to read through the paper, props to you.

  1. The scalar product $(f,g)$ is defined in equation (2.2): $$ (f,g)=\int_{0^-}^{\infty} w(\lambda) \, f(\lambda) \, g(\lambda) \,d\lambda . $$ Note that it includes the weight function $w(\lambda)$.

  2. The polynomials $p_0(\lambda), p_1(\lambda), p_2(\lambda),p_3(\lambda),\dots$ are what you obtain if you start with the monomials $1,\lambda,\lambda^2,\lambda^3,\dots$ and perform the Gram–Schmidt orthogonalization process using this scalar product. This means that they satisfy the orthonormality conditions that $(p_k,p_k)=1$ for each $k$ and $(p_i,p_j)=0$ if $i\neq j$, and moreover each $p_n(\lambda)$ has degree $n$ with positive $\lambda^n$-coefficient.

  3. I'm not 100% sure what you mean by discrete $\lambda$, but if you're referring to the case when the weight function isn't actually a function but a finite linear combination of Dirac delta distributions at the points $\lambda_1,\dots,\lambda_N$, with some weights $w_1,\dots,w_N$, i.e., $$ w(\lambda)=\sum_{m=1}^N w_m \delta_{\lambda_m}(\lambda) , $$ then everything is the same except that the scalar product simplifies to a finite sum instead of an integral, as in equation (2.22): $$ (f,g) = \sum_{m=1}^N w_m \, f(\lambda_m) \, g(\lambda_m) . $$ Moreover, in this case you only get a finite sequence of polynomials $p_0(\lambda),\dots,p_{N-1}(\lambda)$ instead of an infinite sequence. The variable $\lambda$ is still a continuous variable; it's the coefficients of the polynomials $p_n(\lambda)$ that will be functions of the constants $\{ \lambda_m, w_m \}_{m=1}^N$. (See equation (2.38) for an explicit formula.)

  • $\begingroup$ Thank you for this very informative answer. I think I'm understanding it better now. Another further question: In Section 2.2 of the paper it initializes the algorithm with $p_0(\lambda)=1/\sqrt{s_0}$. Since $s_0 = \sum_{m=0}^N w_m$, then is $p_0$ a constant for all $\lambda$? $\endgroup$ – Darcy Apr 20 '16 at 17:26
  • 1
    $\begingroup$ Yes, $p_n$ is of degree $n$ for each $n$, so in particular $p_0$ is of degree zero, i.e., a constant. $\endgroup$ – Hans Lundmark Apr 20 '16 at 17:27
  • $\begingroup$ And finally: how are the subscripts $k$ and $n$ related? Why does he use different subscripts? It seems like they refer to the same index... $\endgroup$ – Darcy Apr 20 '16 at 17:30
  • 1
    $\begingroup$ Looking quickly, I can't see where the subscript $k$ is used. But you can of course use whatever letter you want for the subscript if you want to denote "an arbitrary polynomial in the sequence $p_0,p_1,p_2,\dots$"; call it $p_k$, $p_n$, $p_i$, $p_j$, or anything else, as long as the subscript isn't occupied for something else already. (For example, in the finite case where $N$ is the number of point masses, then $p_{N-1}$ isn't just any of the polynomials, it's specifically the last one.) $\endgroup$ – Hans Lundmark Apr 20 '16 at 17:34

I have had a glance at the paper. It is not that it is badly written but I consider pages 568-573 constitute a very intricated lecture on orthogonal polynomials, and the attached continued fraction. It is not surprising that you cannot easily find your way in this jungle.

The presentation can be vastly simplified by centering all on the positive definite moment matrix and its Cholesky factorization, instead of introducing it at the end.

If you desire, I can send you a very condensed version of my lectures on the subject (20 powerpoint slides or so, in easily understandable French).

Edit Answers to some of your questions:

  • Points i) and iii): parentheses mean dot product in the functional sense $(p_i,p_j)=\int _a^bp_i(\lambda)p_j(\lambda)d \lambda$

  • $n$ is at the same time the index AND the degre of polynomial $p_n$.

  • the discretization made is not classical in my eyes, but I should spend more time on the article.

  • $\begingroup$ Hi, thanks for the response. I am a geophysicist and not a mathematician, so I am not so familiar with a lot of the mathematical background to this paper. I unfortunately don't speak French, but you could send the slides over anyway and I will have a look. My email is dcordell@ualberta.ca $\endgroup$ – Darcy Apr 19 '16 at 20:10
  • $\begingroup$ Also, if you know of any good beginner textbooks on the subject, that would be helpful too. $\endgroup$ – Darcy Apr 19 '16 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.