text:
\section{Introduction}
Numerical integration can be formulated as a finite matrix approximation
of the multiplication operator
\cite{Gragg:1993,
  Simon:2009,
  Velazquez:2008,
  Bultheel+Cantero+Cruz-Barroso:2015,
  Sarmavuori+Sarkka:2019}.
In this article, as in our previous article \cite{Sarmavuori+Sarkka:2019} on this subject, we consider
numerical integration of a composition of two real functions $f(g(\mat{x}))$
as a function of a finite matrix $f(\mat{M}_n[g])$:
\begin{equation}
  \label{eq:numint-def}
  \int_{\Omega} f(g(\mat{x}))\,w(\mat{x})\,d\mat{x}=
  \langle 1,f(\operator{M}[g])\,1\rangle\approx
  \left[
    f(\mat{M}_n[g])
    \right]_{0,0},
\end{equation}
where $\Omega\subset\mathbb{R}^d$, $w(\mat{x})\geq 0$,
$%
  \int_\Omega w(\mat{x})\,d\mat{x}=1
$, %
  $\langle \cdot,\cdot \rangle$ is an inner product
$%
  \langle \phi,\psi \rangle =
  \int_{\Omega}\overline{\phi(\mat{x})}\,\psi(\mat{x})\,
  w(\mat{x})\,d\mat{x}
  $ which defines a norm
  $\|\phi\|=\sqrt{\langle \phi, \phi \rangle}$
  so that
  $\vectorspace{L}^2_w(\Omega)=\{\phi : \|\phi\|<\infty\}$
  is a complete Hilbert space, %
 $\operator{M}[g]$ is a multiplication operator of multiplying with
  a real function $g$, that is, $\operator{M}[g]\,\phi=g\,\phi$
  almost everywhere
  for all
$\phi\in\vectorspace{D}(\operator{M}[g])=\{\phi : \|g\,\phi\|<\infty\}$,
and the elements of the finite matrix approximation
$\mat{M}_n[g]$ of the multiplication operator $\operator{M}[g]$ are
computed in finite dimensional subspace, which is spanned by orthonormal
functions $1=\phi_0,\phi_1,\phi_2,\ldots,\phi_n$, as
\begin{equation}
  \label{eq:matrix-elements}
  \left[
    \mat{M}_n[g]
    \right]_{i,j}
  =\langle \phi_i,\operator{M}[g]\,\phi_j\rangle=\langle \phi_i,g\,\phi_j\rangle
  =
  \int_\Omega \overline{\phi_i(\mat{x})}\,g(\mat{x})\,\phi_j(\mat{x})\,
w(\mat{x})\,d\mat{x}.
\end{equation}

More generally, we can use the same approach to approximate integrals
of the form
\begin{equation*}
  \int_{\Omega} \overline{\phi(\mat{x})}\,f(g(\mat{x}))\,\psi(\mat{x})\,
  w(\mat{x})\,d\mat{x}
  \approx
  \mat{v}_n^*\,f(\mat{M}_n[g])\,\mat{u}_n,
\end{equation*}
where, as $n\rightarrow\infty$,
$\mat{v}_n\rightarrow \mat{v}$ and $\mat{u}_n\rightarrow \mat{u}$
such that
$%
[\mat{v}]_i
=
\langle \phi_i, \phi\rangle$ and %
$        [\mat{u}]_i
        %
        = \langle \phi_i, \psi\rangle
$ %
or integrals
that
involve products
of functions, that is,
\begin{equation*}
    \int_{\Omega} f_1(g_1(\mat{x}))\,f_2(g_2(\mat{x}))\,w(\mat{x})\,d\mat{x}
    \approx
    [f_1(\mat{M}_n[g_1])\,f_2(\mat{M}_n[g_2])]_{0,0}
\end{equation*}
or products of several functions
\begin{equation*}
    \int_{\Omega} \prod_{i=0}^m f_i(g_i(\mat{x}))\,w(\mat{x})\,d\mat{x}
    \approx
    \left[
      \prod_{i=0}^m f_i(\mat{M}_n[g_i])
    \right]_{0,0}.
\end{equation*}

In \cite{Sarmavuori+Sarkka:2019}, this method was shown to be a generalization
of the Gaussian quadrature, that is, we get the Gaussian quadrature when
the inside function is $g(x)=x$ and the basis functions
are polynomials. In \cite{Sarmavuori+Sarkka:2019},
the convergence was considered only
for
bounded outside function $f$
and a bounded inside function $g$.
In this article, our goal is to extend the convergence results to unbounded
functions. This allows us to consider improper integrals of functions
that attain infinity at either or both endpoints of the integration interval.
Our primary interest is in integrals, where infinity
is attained at an endpoint that is infinite or at both endpoints $\pm\infty$.
Proof of convergence for this type of integrals is often
considerably more complicated than for bounded functions on bounded
intervals.
See
\cite{Uspensky:1928,Jouravsky:1928,Shohat+Tamarkin:1963,Freud:1971,
  Bultheel+Diaz-Mendoza+Gonzalez-Vera+Orive:1997,
  Bultheel+Diaz-Mendoza+Gonzalez-Vera+Orive:2000}
for examples of proofs of convergence on unbounded intervals
for Gaussian quadratures
and Gaussian quadratures generalized on rational functions.
Similar results have been proved also purely operator theoretically
for Gaussian quadratures on unbounded intervals for polynomially bounded
functions \cite[Sections~1.8~and~1.13]{Simon:2009}.
We generalize this operator theoretic approach
on larger family of inside functions $g$
and basis functions
and under some conditions
to faster than polynomially growing outside functions $f$.
We also consider improper integrals, where infinity is
attained at a finite endpoint.

As a starting point we take the strong resolvent convergence
(or generalized strong convergence)
\cite{Kato:1995,Reed+Simon:1981,Segal+Kunze:1978,Weidmann:1980,Simon:2015}
that gives well-known convergence results for continuous functions and
characteristic functions.
We extend these results to Riemann--Stieltjes integrable functions
that are by definition bounded but can have infinitely many discontinuities.
We also prove convergence results for a class of linearly bounded functions,
that is, functions that are less than $a+b\,\abs{x}$ for some $a,b>0$.
These results apply to general self-adjoint operators, not just
multiplication operators and therefore they also apply for more exotic
integrals that can be formulated as function of self-adjoint operator
\cite{Segal:1965}.

For Gaussian quadratures, the convergence proofs are often based on
an inequality of the following type (e.g.
\cite{Uspensky:1928,Jouravsky:1928,
Bultheel+Diaz-Mendoza+Gonzalez-Vera+Orive:2000},
\cite[Chapter~IV]{Shohat+Tamarkin:1963} or
\cite[Chapter~3, Lemma~1.5]{Freud:1971}):
\begin{equation}
  \sum_{k=0}^n w_k\,x_k^{2\,m}\leq \int_{-\infty}^\infty x^{2\,m}\,w(x)\,dx,
\label{eq:wx_ineq}
\end{equation}
that is, a Gaussian quadrature approximation with weights $w_k$ and nodes $x_k$
never exceeds the true value of the integral for even monomials $x^{2\, m}$.
We do not have this type of inequality in general but in special cases when
all matrix elements are positive. We give examples of such cases.
In those cases the convergence can still be proved for certain fast growing
classes of functions such as exponentially bounded functions.

A similar inequality to \eqref{eq:wx_ineq} also arises from Jensen's operator inequality when the
outside function $f$ is an operator convex function.
We use this type of inequality for proving convergence of
improper integrals, where infinity of the integrand is attained at a
finite endpoint of the integration interval.

The main contributions of this article are:
\begin{enumerate}
\item
  We show that the theory of strong resolvent convergence can be used
  for analysing the convergence of matrix method of numerical integration
  when the inside function is unbounded.
\item
  We extend the theory of convergence for functions of self-adjoint operators
  in a converging sequence from continuous and bounded functions to
  wider classes of discontinuous and unbounded functions.
\item
  We show theoretically and by a numerical example that convergence holds
  for wider class of unbounded functions when the matrix elements are
  nonnegative.
\item
  We also show the convergence of improper integrals on a finite endpoint for
  functions bounded by an operator convex function by using Jensen's operator
  inequality.
\end{enumerate}

\section{Preliminaries}
\label{sec:prelim}

A self-adjoint operator is such operator that $\operator{A}=\operator{A}^*$.
When a function $g$ is a measurable real function,
the maximal multiplication operator
$\operator{M}[g]$, that is,
the multiplication operator with domain
$\vectorspace{D}(\operator{M}[g])=\{\phi : \|g\,\phi\|<\infty\}$,
is self-adjoint.

We denote strong convergence as
$\operator{A}_n\xrightarrow{s}\operator{A}$. Strong resolvent
convergence of a sequence of self-adjoint operators $\operator{A}_n$
to a self-adjoint operator $\operator{A}$ means that
$(\operator{A}_n - z)^{-1}\xrightarrow{s}(\operator{A}-z)^{-1}$
for some nonreal $z$
which
we denote as $\operator{A}_n\xrightarrow{srs}\operator{A}$
\cite{Kato:1995,Reed+Simon:1981,Segal+Kunze:1978,Weidmann:1980,Simon:2015}.
For a function
of a self-adjoint operator, we have the following useful theorem
about the convergence of functions of operators.
\begin{theorem}
  \label{thm:bounded-function-srs}
  If a sequence of self-adjoint operators $\operator{A}_n$ converges in strong
  resolvent sense to a self-adjoint operator $\operator{A}$, and $f$ is a bounded
  continuous function defined on $\mathbb{R}$, then $f(\operator{A}_n)$
  converges strongly to $f(\operator{A})$.
\end{theorem}
\begin{proof}
  See
  \cite[Theorem~VIII.20 (b)]{Reed+Simon:1981},
  \cite[Theorem~9.17]{Weidmann:1980}
  or \cite[Theorem~11.4]{Segal+Kunze:1978}.
\end{proof}

The spectral representation of a self-adjoint operator is
$%
  \operator{A} = \int_{\sigma(\operator{A})}t\,d\operator{E}(t)
$ %
where $\sigma(\operator{A})$ is the spectrum, and $\operator{E}(t)$ is the
spectral family of the operator $\operator{A}$
\cite{Stone:1932, Riesz+Nagy:1955,Akhiezer+Glazman:1993,Kato:1995,
  Halmos:1982, Kreyszig:1989, Reed+Simon:1981, Segal+Kunze:1978,
  Weidmann:1980, Schmudgen:2012, Simon:2015}.
For any function which is measurable
with respect to the spectral family, we have
\begin{equation}
  \label{eq:spectral-rep}
  f(\operator{A})=\int_{\sigma(\operator{A})}f(t)\,d\operator{E}(t).
\end{equation}
The spectral family can be defined as left or right
continuous with respect to the parameter $t$ in strong convergence sense
\cite[Chapter~VI, Section~5.1]{Kato:1995}.
We adopt the usual convention of defining it as right continuous. Only points
of discontinuity of $\operator{E}(t)$ are the eigenvalues of $\operator{A}$
\cite[Theorem~7.23]{Weidmann:1980}. An \textit{eigenvalue}
of an operator is a value $\lambda\in\mathbb{C}$ that satisfies
$%
  \lambda\,\phi=\operator{A}\,\phi.
$ %
  Function $\phi$ is called an \textit{eigenfunction} or an eigenvector.
  We denote the set of eigenvalues of operator $\operator{A}$ as
  $\sigma_p(\operator{A})$.

An eigenvalue $\lambda$ of a multiplication operator $\operator{M}[g]$ has to
satisfy
\cite[p.~103, Example~1]{Weidmann:1980}
\begin{equation*}
  \int_{\{\mat{x}\in\Omega\,:\,\lambda=g(\mat{x})\}}w(\mat{x})\,d\mat{x}>0
\end{equation*}
and an eigenfunction has to be scalar multiple of a characteristic function
$\phi(\mat{x})=c\,\chi_S(\mat{x})$
for $c\in\mathbb{C}$ and some set
$S\subset\{\mat{x}\in\Omega : \lambda=g(\mat{x})\}$,
that is, measurable with respect to the weight function $w(\mat{x})$.
The spectrum of the multiplication operator $\sigma(\operator{M}[g])$
for a real function $g$
is the essential range of the function $g$
\cite[Problem~67]{Halmos:1982},
\cite[Chapter~VIII.3, Proposition~1]{Reed+Simon:1981},
\cite[p.~103, Example~1]{Weidmann:1980},
that is,
\begin{equation*}
  \sigma(\operator{M}[g])=
  \mathcal{R}(g)=
  \left\{
  y\in\mathbb{R} :
  \text{for all }\epsilon>0,
  \int_{\{\mat{x}\in\Omega\,:\,|g(\mat{x})-y|<\epsilon\}}
  w(\mat{x})\,d\mat{x}>0
  \right\}.
\end{equation*}

For the convergence of the spectral families, we have the following theorem.
\begin{theorem}
  \label{thm:characteristic-convergence}
  Let $\operator{A}_n$ and $\operator{A}$ be self-adjoint operators with
  spectral families $\operator{E}_n(t)$ and $\operator{E}(t)$, respectively.
  If $\operator{A}_n\xrightarrow{srs}\operator{A}$, then we have the following
  equivalent results:
  \begin{enumerate}
    \item $\operator{E}_n(t)\xrightarrow{s}\operator{E}(t)$ when $t$ is
      not an eigenvalue of $\operator{A}$.
    \item For characteristic functions we have
      $\chi_{[a,b]}(\operator{A}_n)\xrightarrow{s}\chi_{[a,b]}(\operator{A})$
      when $a$ and $b$ are not eigenvalues of $\operator{A}$.
  \end{enumerate}
\end{theorem}
\begin{proof}
  For 1, see \cite[Chapter~VIII, Theorem~1.15]{Kato:1995} or
  \cite[Theorem~9.19]{Weidmann:1980}. For 2, see
  \cite[Theorem~VIII.24 (b)]{Reed+Simon:1981},
  \cite[Theorem~11.4 e]{Segal+Kunze:1978}, \cite[Theorem~9.19]{Weidmann:1980}
  or \cite[Chapter~7.2, problem~5]{Simon:2015}.
\end{proof}
This theorem has two important consequences: First, it shows that there is
a converging measure and, in fact, several of them, not just
$\langle 1,\operator{E}_n(t)\,1\rangle
\rightarrow \langle 1,\operator{E}(t)\,1\rangle$, but all
$\langle \psi,\operator{E}_n(t)\,\phi\rangle
\rightarrow \langle \psi,\operator{E}(t)\,\phi\rangle$.
Second, it shows that we cannot expect convergence for all functions
in \eqref{eq:spectral-rep} that are measurable with respect to
$\operator{E}(t)$. The convergence can fail if the function $f$
 is discontinuous at an eigenvalue of $\operator{A}$.

A self-adjoint operator has an infinite matrix representation with
 matrix elements \eqref{eq:matrix-elements} if and only if
 the vectors
of the orthonormal
basis
$\phi_0,\phi_1,\ldots$ are dense in the underlying Hilbert space, and
\begin{equation}
  \label{eq:infinite-matrix-condition}
  \|\operator{A}\,\phi_i\|<\infty
\end{equation}
for all $i=0,1,\ldots$ (see \cite[Theorem~3.4]{Stone:1932} or
\cite[Section~47]{Akhiezer+Glazman:1993}).
For a multiplication operator $\operator{M}[g]$ this
is equivalent to
\begin{equation}
  \label{eq:mult-mat}
  \int_{\Omega} |g(\mat{x})\,
  \phi_i(\mat{x})|^2\,w(\mat{x})\,d\mat{x}<\infty
\end{equation}
for all $i=0,1,\ldots$.
If $\operator{A}$ is a self-adjoint operator with an infinite matrix
representation $\mat{A}_\infty$, then also $\mat{A}_\infty$ is a self-adjoint
operator on $\ell^2$, the space of absolutely square-summable sequences.
The domain of $\mat{A}_\infty$ is
$
\vectorspace{D}(\mat{A}_\infty)=
  \left\{
  \mat{v} : \|\mat{A}_\infty\,\mat{v}\|<\infty
  \right\}
$
  %
which is isomorphic to $\vectorspace{D}(\operator{A})$.

\subsection{Matrices with nonnegative elements}
\label{sec:matrices-with-nonnegative-elements}
We define an elementwise partial order of matrices as
  $\mat{A}{\scriptstyle\overset{ew}\leq}\mat{B}$,
  that is, $\mat{A}{\scriptstyle\overset{ew}\leq}\mat{B}$ if $[\mat{A}]_{i,j}\leq [\mat{B}]_{i,j}$ for all
  $i,j$.
For functions of finite matrices with nonnegative elements, we introduce
two classes of special interest \cite{Hansen:1992,Hansen:2000}.
\begin{definition}
  \label{def:m-positive}
  Let $f:I\mapsto\mathbb{R}$ be a real function defined on an interval
  $I\subset\mathbb{R}$.
  \begin{enumerate}
  \item
    $f$ is \textit{$m$-positive} if $0\in I$ and $f(\mat{A}){\scriptstyle\overset{ew}\geq} \mat{0}$ for every
    symmetric $n\times n$ matrix $\mat{A}{\scriptstyle\overset{ew}\geq} \mat{0}$ with spectrum contained
    in $I$, and every $n\in \mathbb{N}$.
  \item
    $f$ is \textit{$m$-monotone} if
    $\mat{A}{\scriptstyle\overset{ew}\leq} \mat{B}\Rightarrow f(\mat{A}){\scriptstyle\overset{ew}\leq} f(\mat{B})$,
    for all symmetric $n\times n$ matrices $\mat{A},\mat{B}{\scriptstyle\overset{ew}\geq} \mat{0}$
      with spectra in $I$, and every $n\in\mathbb{N}$.
  \end{enumerate}
\end{definition}
Under favorable conditions, the concepts of $m$-positive and $m$-monotone
functions are equal as is
shown by the following theorem.
\begin{theorem}
  \label{thm:m-positive-equal-m-monotone}
  Let $f:I\mapsto\mathbb{R}$ be a real function defined on an interval
  $I\subset\mathbb{R}$ such that $0\in I$ and $f(0)\geq 0$. Then
  $f$ is $m$-positive if and only if it is $m$-monotone.
\end{theorem}
\begin{proof}
  See \cite[Theorem~2.1~(1)]{Hansen:1992}.
\end{proof}
The usefulness of the theorem stems from the fact that it is often easier
to see that a function is $m$-positive than to directly see that it is
$m$-monotone.
On $[0,\infty)$ and $(-\infty,\infty)$ intervals,
all $m$-positive and hence $m$-monotone functions are absolutely
monotone functions, that is, functions that have a converging
power series expansion $f(x)=\sum_{n=0}^\infty c_n\,x^n$ for
$n=0,1,2,\ldots$ and $c_n\geq 0$ for all $n$
\cite{Micchelli+Willoughby:1979,Hansen:1992,Hansen:2000}.

\subsection{Jensen's operator inequality}
\label{sec:jensens-operator-inequality}
In this section, we use a different concept of nonnegativity than in the
previous section which hence also imposes a different partial order.
Instead of elementwise nonnegativity,
for bounded operators, we say that a self-adjoint operator
$\operator{A}$ is nonnegative
if
$
  0 \leq \langle \phi, \operator{A}\,\phi\rangle
$
  %
for all
$\phi$ in the Hilbert space. We mark
this as $\operator{A}{\scriptstyle\overset{qf}\geq} 0$.
For bounded operators, we can use definition $\operator{A}{\scriptstyle\overset{qf}\leq} \operator{B}$
if $0 {\scriptstyle\overset{qf}\leq} \operator{B}-\operator{A}$
\cite{Kato:1995,Reed+Simon:1981,Bhatia:1997,Carlen:2010}.

We notice that for bounded self-adjoint operators
$\operator{A}{\scriptstyle\overset{qf}\geq} 0 \iff \sigma(\operator{A})\subset [0,\infty)$.
  Therefore, we generalize the nonnegativity to unbounded self-adjoint
  operators so that
  an unbounded self-adjoint operator $\operator{A}{\scriptstyle\overset{qf}\geq} 0$ if
  $\sigma(\operator{A})\subset [0,\infty)$.
    %
    %
    %
    %
    %
    %
    We want to extend the order relation of self-adjoint operators
    to unbounded operators naturally, so
    that for self-adjoint operators $\operator{A},\operator{B}$, we have
    $\operator{A}{\scriptstyle\overset{qf}\leq}\operator{B}$ if
    $0{\scriptstyle\overset{qf}\leq}\operator{B}-\operator{A}$. The definition of operator
    $\operator{B}-\operator{A}$ might be complicated in general. For
    unbounded self-adjoint infinite matrices $\mat{A}_\infty,\mat{B}_\infty$,
    the difference is possible to simply define elementwise as
    $[\mat{B}_\infty-\mat{A}_\infty]_{i,j}=
    [\mat{B}_\infty]_{i,j}-[\mat{A}_\infty]_{i,j}$.
    The elementwise difference of infinite matrices satisfies
    \eqref{eq:infinite-matrix-condition}, and the domain is simply
    $\vectorspace{D}(\mat{B}_\infty-\mat{A}_\infty)=
    \{\mat{v} : \|(\mat{B}_\infty-\mat{A}_\infty)\,\mat{v}\|<\infty\}$.
    Then we can say that $\mat{B}_\infty{\scriptstyle\overset{qf}\geq}\mat{A}_\infty$
    if $\mat{B}_\infty-\mat{A}_\infty{\scriptstyle\overset{qf}\geq} 0$.

The operator convexity is usually defined for bounded operators
\cite{Bhatia:1997,Hansen+Pedersen:2003,Carlen:2010}.
\begin{definition}
  A real function $f$ defined on an interval $I$ is said to be
  operator convex for bounded operators
  if for all self-adjoint bounded operators
  $\operator{A},\operator{B}$ %
  with spectrum in $I$,
  for each $\lambda\in[0,1]$, we have
  \begin{equation*}
    f(\lambda\,\operator{A}+(1-\lambda)\,\operator{B}){\scriptstyle\overset{qf}\leq} \lambda\,
    f(\operator{A})+(1-\lambda)\,f(\operator{B}).
  \end{equation*}
\end{definition}
In the above definition, also $f(\operator{A})$ and $f(\operator{B})$
are bounded operators because operator convex functions are continuous.

The set of operator convex functions does not contain all convex functions.
For example, for bounded operators on interval
$I=(0,\infty)$, $f(x)=x^p$ is operator convex for
  $-1\leq p\leq 0$ and $1\leq p \leq 2$ but not for
  $p<-1$, $0 < p < 1$ or $2<p$
\cite[Theorem~2.6]{Carlen:2010}. The set of operator convex functions
can be even smaller for unbounded operators than it is for bounded operators.
On nonnegative interval $[0,\infty)$, it has been shown that operator convex
  functions for
bounded operators are also operator convex for unbounded operators
\cite{Dinh+Tikhonov+Veselova:2019}.

For the proof of convergence of the matrix approximation for some integrals,
we can use the following Jensen's operator inequality which we will
extend for unbounded self-adjoint infinite matrices.
\begin{theorem}
  \label{thm:Jensen}
  For a continuous and operator convex function $f$ defined on an interval
  $I$
  \begin{equation*}
    \operator{P}\,
    f(\operator{P}\,\operator{A}\,\operator{P}+s\,(1-\operator{P}))\,
    \operator{P}{\scriptstyle\overset{qf}\leq} \operator{P}\,f(\operator{A})\,\operator{P}
  \end{equation*}
  for each projection operator $\operator{P}$ on an infinite dimensional
  Hilbert space, every bounded
  self-adjoint operator $\operator{A}$ with spectrum
  in $I$ and every $s\in I$.
\end{theorem}
\begin{proof}
  See \cite[Theorem~2.1 (i) and (iv)]{Hansen+Pedersen:2003}.
\end{proof}

We extend the definition of operator convexity
to nonnegative functions of unbounded self-adjoint infinite matrices.
A nonnegative function of an operator is
always a nonnegative operator.
For each nonnegative operator,
there exists
a corresponding element in the extended positive cone of the von Neumann
algebra of the
bounded
operators on the underlying Hilbert space.
The convex combination of elements of the extended positive cone
is well defined. Thus, we can define
operator convex functions so that the domain of the function
is unbounded operators with infinite matrix representation and
the nonnegative function of an operator in that domain is an element in the
extended positive cone.

We present the definitions of von Neumann algebra and extended positive
cone of von Neumann algebra as in \cite{Takesaki:2001,Takesaki:2002}.
\begin{definition}
Commutant
$\mathcal{M}^\prime$ of
a subset $\mathcal{M}$ of the bounded operators is a set of all bounded
operators that commute with every operator in $\mathcal{M}$
\cite[Chapter~2, Section~3]{Takesaki:2001}.
A von Neumann algebra
$\mathcal{M}$ is such a subalgebra of bounded operators on a Hilbert space
that $\mathcal{M}$ is a commutant of its own commutant, that is,
$\mathcal{M}=\mathcal{M}^{\prime\prime}$ \cite[Theorem~II.3.9]{Takesaki:2001}.
The space $\mathcal{M}_*$
is the Banach space of all $\sigma$-weakly continuous linear functionals
on von Neumann algebra $\mathcal{M}$, also called normal linear functionals
\cite[Theorem~III.2.7]{Takesaki:2001}.
Space $\mathcal{M}_*^+$ is the space of all positive normal linear functionals,
that is, such functionals $\varphi$ that
$\varphi(\operator{A}^*\,\operator{A})\geq 0$
for all $\operator{A}\in\mathcal{M}$ \cite[Definition~I.9.4]{Takesaki:2001}.
The extended positive cone $\widehat{\mathcal{M}}_+$ of a von Neumann algebra
$\mathcal{M}$ is the set of maps $m:\mathcal{M}_*^+\mapsto[0,\infty]$ that
satisfy (\cite[Definition~IX.4.4]{Takesaki:2002} or
\cite[Definition~2.2]{Dinh+Tikhonov+Veselova:2019})
\begin{enumerate}
\item $m(\lambda\,\varphi)=\lambda\,m(\varphi),~\lambda\geq 0$
  (here $0\cdot(+\infty)=0$),
\item $m(\varphi+\psi)=m(\varphi)+m(\psi)$ for $\varphi,\psi\in\mathcal{M}_*^+$,
\item $m$ is lower semicontinuous.
\end{enumerate}
\end{definition}
For example, let $\mathcal{M}$ be the von Neumann algebra of all bounded
operators
on some Hilbert space, unbounded self-adjoint operator $\operator{A}{\scriptstyle\overset{qf}\geq} 0$
and vector
$\phi\in\vectorspace{D}(\operator{A})$. Then
there are $\varphi_\phi\in\mathcal{M}_*^+$ and
$m_{\operator{A}}\in\widehat{\mathcal{M}}_+$ such that
$\langle\phi,\operator{A}\,\phi\rangle=\varphi_\phi(\operator{A})
=m_{\operator{A}}(\varphi_\phi)$.
We can define order relation of two unbounded self-adjoint operators
$\operator{A},\operator{B}{\scriptstyle\overset{qf}\geq} 0$ so that $\operator{A}{\scriptstyle\overset{qf}\geq}\operator{B}$
if the corresponding $m_{\operator{A}},m_{\operator{B}}\in\widehat{\mathcal{M}}_+$
satisfy $m_{\operator{A}}(\varphi)\geq m_{\operator{B}}(\varphi)$ for all
$\varphi\in\mathcal{M}_*^+$.

\section{Convergence}
\label{sec:convergence}
We establish convergence for growing class of outside functions $f$. We start
by first proving the strong resolvent convergence that then immediately
covers the
bounded continuous outside functions in
Section~\ref{sec:proof-of-strong-resolvent-convergence}.
Our first extension to this basic result is the class of discontinuous
functions in Section~\ref{sec:convergence-for-discontinuous-functions}.

Then we extend the results to unbounded functions.
First, we discuss the topic of general unbounded outside functions in
Section~\ref{sec:strong-resolvent-convergence-for-unbounded-functions}
and notice that it is too wide class of functions to prove convergence
in general. Then we find proof for growing class of unbounded functions
as the inside function $g$ and basis functions $\phi_i$ are limited.
In Section~\ref{sec:convergence-for-quadratically-bounded-functions},
the convergence is proved for quadratically bounded outside functions
without restrictions to the inside function and polynomially bounded
outside functions when the inside function and the basis functions are
polynomials.
In Section~\ref{sec:nonnegative-matrix-coefficients},
the convergence is proved for inside functions that have matrix
representations with nonnegative coefficients. The outside functions do not
even have to meet the criteria of being functions of operators that have
infinite matrix representation.

Finally, we consider convergence for outside functions that are singular on a
finite endpoint of an integration interval in
Section~\ref{sec:improper-integrals-on-finite-endpoints-of-an-interval}.

\subsection{Proof of strong resolvent convergence}
\label{sec:proof-of-strong-resolvent-convergence}
It can be difficult to prove directly that a sequence of self-adjoint
operators converges in strong resolvent sense. Therefore, we will
introduce a concept of core
that can be used for proving the strong
resolvent convergence.
\begin{definition}
A core $\vectorspace{D}_0$ of
an operator $\operator{A}$ is a subspace of the domain of the operator
such that the second adjoint of the restriction of the operator to the core is
the
operator itself, that is,
$
(\operator{A}|_{\vectorspace{D}_0})^{**}=\operator{A}.
$
\end{definition}
Usually the core is defined in terms of the closure of the operator
\cite{Kato:1995, Reed+Simon:1981, Weidmann:1980, Simon:2015}.
In this article, we note that for densely defined and closable operators,
the closure is equivalent to the second adjoint
\cite[Chapter~III, Theorem~5.29]{Kato:1995},
\cite[Theorem~VIII.1~(b)]{Reed+Simon:1981},
\cite[Theorem~5.3~(b)]{Weidmann:1980}
or
      \cite[Theorem~7.1.1~(c)]{Simon:2015}.

With a suitable core, the strong resolvent convergence can be proved by
the following theorem.
\begin{theorem}
  Let $\vectorspace{D}_0$ be a core of self-adjoint operators $\operator{A}_n$
  and $\operator{A}$ where $n = 0,1,2,\ldots$. If
  for each $\phi\in\vectorspace{D}_0$,
  \begin{equation*}
  \lim_{n\rightarrow \infty}\|(\operator{A}-\operator{A}_n)\,\phi\|=0,
  \end{equation*}
  then $\operator{A}_n\xrightarrow{srs}\operator{A}$.
\end{theorem}
\begin{proof}
  See   \cite[Theorem~VIII.25~(a)]{Reed+Simon:1981} for proof, or
  \cite[Chapter~VIII, Corollary~1.6]{Kato:1995},
  \cite[Theorem~9.16~(i)]{Weidmann:1980}, or
  \cite[Theorem~7.2.11]{Simon:2015} for proofs of similar but slightly more
  general theorems.
\end{proof}
After finding a suitable core, we can prove the strong resolvent convergence. The space of absolutely square summable sequences or infinite dimensional vectors is $\ell^2$. Its subspace $\ell^2_0$ is a space, where vectors have only finitely many non-zero elements. It now turns out that $\ell^2_0$ is in general a core for any self-adjoint infinite matrix.
Although similar ideas related to infinite band matrices have been discussed in
\cite{Barrios+Lopez+Martinez-Finkelshtein+Torrano:1999,Dombrowski:1990} and
a space which is isomorphic to $\ell^2_0$ has been used as a core in
\cite{Arai:1992,Hansen:2008}, it looks like
a proof has been presented for a space that is isomorphic
to $\ell^2_0$ only in a special case
\cite[Lemma~4.4]{Rosler:2019}, so
we present a proof in the general case.
\begin{theorem}
  Let $\operator{A}$ be a self-adjoint operator with an infinite
  matrix representation $\mat{A}_\infty$.
  Then $\ell^2_0$ is a core for $\mat{A}_\infty$.
\end{theorem}
\begin{proof}
  Our proof is an adoption of proof of \cite[Theorem~6.20]{Weidmann:1980}.
  We use the following property of the adjoint:
  \begin{itemize}
    \item Let $\operator{A}$ and $\operator{B}$ be densely defined operators,
      then
      $\vectorspace{D}(\operator{A})\subset\vectorspace{D}(\operator{B})$
      and $\operator{A}=\operator{B}|_{\vectorspace{D}(\operator{A})}$
      $\Rightarrow$
      $\vectorspace{D}(\operator{B}^*)\subset\vectorspace{D}(\operator{A}^*)$
      and $\operator{B}^*=\operator{A}^*|_{\vectorspace{D}(\operator{B}^*)}$,
      that is,
      \begin{equation}
        \label{eq:operator-subset}
        \operator{A}\subset\operator{B}
        \Rightarrow
        \operator{B}^*\subset\operator{A}^*.
      \end{equation}
      (see \cite[p.~72]{Weidmann:1980} or \cite[p.~252]{Reed+Simon:1981}.)
  \end{itemize}
  Let operator $\operator{A}_0=\mat{A}_\infty|_{\ell^2_0}$.
  The natural basis vectors of $\ell^2$ are vectors $\mat{e}_i$
  that have 1 in component $i$ and 0 in other components.
  Vectors of $\ell^2_0$ are finite linear combinations of basis vectors
  $\mat{e}_i$.
  Because all basis vectors $\mat{e}_i\in\vectorspace{D}(\mat{A}_\infty)$,
  we also have $\ell^2_0\subset\vectorspace{D}(\mat{A}_\infty)$.
  Therefore, we have $\operator{A}_0\subset\mat{A}_\infty$.
  By \eqref{eq:operator-subset}, we have
  $\mat{A}_\infty^*\subset \operator{A}_0^*$. Because $\mat{A}_\infty$
  is self-adjoint, we have
  $\mat{A}_\infty=\mat{A}_\infty^*\subset\operator{A}_0^*$.
  We show that also
  $\operator{A}_0^*\subset\mat{A}_\infty$, that is,
  %
  $
    \vectorspace{D}(\operator{A}_0^*)
    \subset
    \vectorspace{D}(\mat{A}_\infty)
    =
    \{\mat{u} : \|\mat{A}_\infty\,\mat{u}\|^2<\infty\}
    $
    %
  and
  $\operator{A}_0^*=\mat{A}_\infty|_{\vectorspace{D}(\operator{A}_0^*)}$.
  Let $\mat{u}\in\vectorspace{D}(\operator{A}_0^*)$ and thus
  $\|\operator{A}_0^*\,\mat{u}\|^2<\infty$.
  Because $\mat{e}_k\in\vectorspace{D}(\operator{A}_0)$, we have for all
  $k=0,1,\ldots$
  \begin{equation}
    \label{eq:el-square-zero-dual}
    \langle \mat{e}_k,\operator{A}_0^*\,\mat{u}\rangle
    =
    \langle \operator{A}_0\,\mat{e}_k,\mat{u}\rangle
    =\sum_{j=0}^\infty \overline{a_{j,k}}\,u_j
    =\sum_{j=0}^\infty a_{k,j}\,u_j
    =\langle \mat{e}_k,\mat{A}_\infty\,\mat{u} \rangle.
  \end{equation}
  We square and sum over $k$
  \begin{equation*}
    \|\mat{A}_\infty\,\mat{u}\|^2
    =
    \sum_{k=0}^\infty|\langle \mat{e}_k,\mat{A}_\infty\,\mat{u}\rangle|^2
    =
    \sum_{k=0}^\infty|\langle \mat{e}_k,\operator{A}_0^*\,\mat{u}\rangle|^2
    =
    \|\operator{A}_0^*\,\mat{u}\|^2<\infty,
  \end{equation*}
  that is,
  $\mat{u}\in\vectorspace{D}(\operator{A}_0^*)\Rightarrow
  \mat{u}\in\vectorspace{D}(\mat{A}_\infty)$.
  By \eqref{eq:el-square-zero-dual} we also have for any
  $\mat{u}\in\vectorspace{D}(\operator{A}_0^*)$
  \begin{equation*}
    \operator{A}_0^*\,\mat{u}=
    \sum_{k=0}^\infty
    \langle \mat{e}_k,\operator{A}_0^*\,\mat{u}\rangle\,
    \mat{e}_k=
    \sum_{k=0}^\infty
    \langle \mat{e}_k,\mat{A}_\infty\,\mat{u}\rangle\,
    \mat{e}_k=\mat{A}_\infty\,\mat{u},
  \end{equation*}
  that is,
  $\operator{A}_0^*=\mat{A}_\infty|_{\vectorspace{D}(\operator{A}_0^*)}$.
  Because we have $\mat{A}_\infty\subset\operator{A}_0^*$ and
  $\operator{A}_0^*\subset\mat{A}_\infty$, we have
  $\operator{A}_0^*=\mat{A}_\infty$. Therefore, we also have
  $\operator{A}_0^{**}=\mat{A}_\infty^*=\mat{A}_\infty$.
\end{proof}

After discovering the core, it is straightforward to prove the
strong resolvent convergence.
\begin{theorem}
  \label{thm:srs-convergence}
  Let $\operator{A}$ be a self-adjoint operator with an infinite matrix
  representation
$\mat{A}_\infty$. Let finite matrices $\mat{A}_n$ be leading
principal $(n+1)\times (n+1)$
submatrices of $\mat{A}_\infty$. Then
for any $\mat{v}\in\ell^2_0$
\begin{equation*}
  \left[
    \begin{array}{ccc}
    \mat{A}_n & 0 & \ldots\\
    0         & 0 & \ldots\\
    \vdots    & \vdots & \ddots
    \end{array}
    \right]\,\mat{v}
  \rightarrow \mat{A}_\infty\,\mat{v},
\end{equation*}
and
\begin{equation*}
  \left[
    \begin{array}{ccc}
    \mat{A}_n & 0 & \ldots\\
    0         & 0 & \ldots\\
    \vdots    & \vdots & \ddots
    \end{array}
    \right]
  \xrightarrow{srs} \mat{A}_\infty.
\end{equation*}
\end{theorem}
\begin{proof}
  Let $\mat{I}_n$ be $(n+1)\times (n+1)$ identity matrix.
  With the identity matrix, we can write
  \begin{equation*}
    \left[
      \begin{array}{ccc}
        \mat{A}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right]
    =
    \left[
      \begin{array}{ccc}
        \mat{I}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right]
    \, \mat{A}_\infty \,
    \left[
      \begin{array}{ccc}
        \mat{I}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right].
  \end{equation*}
  We take an arbitrary vector in the core $\mat{v}\in\ell^2_0$. Let $m$ be the
  last index, where components of $\mat{v}$ are non-zero, that is,
  $v_i=0$ for all $i>m$. For all $n>m$,
  \begin{equation*}
    \left[
      \begin{array}{ccc}
        \mat{I}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right]
    \, \mat{v}=\mat{v},
  \end{equation*}
  and thus
  \begin{equation*}
    \left(
    \left[
      \begin{array}{ccc}
        \mat{A}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right]
    -\mat{A}_\infty
    \right)\,\mat{v}
    =
    \left(
    \left[
      \begin{array}{ccc}
        \mat{I}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right]
    -\mat{I}_\infty
    \right)
    \,\mat{A}_\infty\,\mat{v}.
  \end{equation*}
  Vector $\mat{u}=\mat{A}_\infty\,\mat{v}\in\ell^2$, and therefore we have
  \begin{equation*}
    \left(
    \left[
      \begin{array}{ccc}
        \mat{I}_n & 0 & \ldots\\
        0         & 0 & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right]
    -\mat{I}_\infty
    \right)
    \,\mat{u}\rightarrow\mat{0}
  \end{equation*}
  as $n\rightarrow\infty$.
\end{proof}

From this we already have convergence for bounded and continuous functions.
\begin{theorem}
  Let orthonormal functions $\phi_0=1,\phi_1,\phi_2,\ldots$ be dense in
  $\vectorspace{L}^2_w(\Omega)$. Also let function
  $g$ satisfy \eqref{eq:mult-mat} and matrices $\mat{M}_n[g]$ have
  elements as in \eqref{eq:matrix-elements}.
  Then for bounded and continuous $f$
\begin{equation}
  \label{eq:bounded-convergence}
  f\left(
  \left[
    \begin{array}{ccc}
      \mat{M}_n[g] & 0 & \ldots\\
      0         & 0 & \ldots\\
      \vdots    & \vdots & \ddots
    \end{array}
    \right]
  \right)
  \xrightarrow{s}
  f(\mat{M}_\infty[g]),
\end{equation}
and
\begin{equation*}
\mat{e}_0^\top\,f(\mat{M}_n[g])\,\mat{e}_0
\rightarrow
\mat{e}_0^\top\,f(\mat{M}_\infty[g])\,\mat{e}_0
=
\langle 1,f(\operator{M}[g]\,1\rangle
=
\int_\Omega f(g(\mat{x}))\,w(\mat{x})\,d\mat{x}.
\end{equation*}
\end{theorem}
\begin{proof}
  By Theorem~\ref{thm:srs-convergence}, we have
  \begin{equation*}
    \left[
      \begin{array}{ccc}
        \mat{M}_n[g] & 0      & \ldots\\
        0            & 0      & \ldots\\
        \vdots       & \vdots & \ddots
      \end{array}
      \right]
        \xrightarrow{srs}\mat{M}_\infty[g],
  \end{equation*}
  and the convergence of \eqref{eq:bounded-convergence} follows from
  Theorem~\ref{thm:bounded-function-srs}.
  The convergence of quadratic form
  $\mat{e}_o^\top\,f(\mat{M}_n[g])\,\mat{e}_0\rightarrow
  \mat{e}_0^\top\,f(\mat{M}_\infty[g])\,\mat{e}_0$ is a special case of
  weak convergence that follows from the strong convergence
  \eqref{eq:bounded-convergence}.
\end{proof}

This result already extends the previous result
\cite[Theorem~3]{Sarmavuori+Sarkka:2019}
from bounded functions $g$ to unbounded functions $g$.
However, our aim is to further extend the results to discontinuous and unbounded functions $f$ as well.

\subsection{Convergence for discontinuous functions}
\label{sec:convergence-for-discontinuous-functions}
First, we consider convergence when the integrand function $f$ in
\eqref{eq:numint-def} is bounded but discontinuous. The largest space of
discontinuous functions that are possible to integrate numerically are the
Riemann--Stieltjes integrable functions.
It does not seem to be
possible to obtain convergence results for larger set of functions
like, for example,
Lebesgue--Stieltjes integrable functions
\cite[Chapter~1.8]{Davis+Rabinowitz:1984}. Similar arguments rule out
convergence for Borel or Baire measurable functions that are not
Riemann--Stieltjes
integrable.
When the starting point is a Borel measure, the largest set of
discontinuous functions which the convergence theorem covers seem to be
 functions
that are discontinuous on a set that is contained on a closed set which has
zero measure with respect to the spectral family of the operator
\cite[Theorem~2.6]{Bade:1954}, \cite{Foguel:1958},
\cite[Theorem~3]{Simpson:1966}
or
\cite[Theorem~3]{Sarmavuori+Sarkka:2019}.
These convergence theorems do not then
cover, for example,
Thomae's function \cite[Example~7.1.7]{Bartle+Sherbert:2011}
which is discontinuous on rational numbers but is
Riemann--Stieltjes integrable.
The definition of the Riemann-- and the related Darboux--Stieltjes integrals
are
 given in \ref{app:integrable-functions}.

The spectral integral is usually defined as a Borel or a Baire measure.
One example of spectral integral as a Riemann--Stieltjes integral
can be found in \cite[Chapter~4.1]{Schmudgen:2012}. However,
the Riemann--Stieltjes integral is defined there in norm operator topology
and we will rather make the definition in strong operator topology because
the finite matrix approximations
can converge in norm only if the operator is compact
\cite[Problem~175]{Halmos:1982},
\cite[Theorem~1.2]{Morrison:1995}
and
the function space defined in the norm operator topology can only be a subset
of the space defined in the strong operator topology.
We define the spectral integral as Darboux--Stieltjes integral
in a same manner as it is defined as Lebesgue--Stieltjes integral
in \cite[Section~126]{Riesz+Nagy:1955}.
The Riemann--Stieltjes integral is then the same as the Darboux--Stieltjes
integral except limited to functions that are continuous at the discontinuities
of the spectral family.
We can define a non-decreasing function $\rho(t)$ through the spectral family
as $\rho(t)=\|\operator{E}(t)\,\phi\|^2$ for any $\phi$ in the Hilbert space.
Because operators $\operator{E}(t)$ are projections,
$\|\operator{E}(t)\,\phi\|^2=\langle \phi,\operator{E}(t)\,\phi\rangle$.
Thus, we can define
\begin{equation*}
  \langle \phi, f(\operator{A})\,\phi\rangle=
  \int_{\sigma(\operator{A})}f(t)\,d\langle \phi,\operator{E}(t)\,\phi\rangle
  =
  \int_{\sigma(\operator{A})}f(t)\,d\|\operator{E}(t)\,\phi\|^2
\end{equation*}
as a Darboux--Stieltjes integral. We say that a function $f$ is
Darboux--Stieltjes
integrable with respect to $\|\operator{E}(t)\,\phi\|^2$ if the
Darboux--Stieltjes integral exists.

We can use polarization identity \cite[Section~126]{Riesz+Nagy:1955}
\begin{equation*}
  \langle \phi, \operator{E}(t)\,\psi\rangle
  =
  \left\|
  \operator{E}(t)\,\frac{\scriptstyle\phi+\psi}{\scriptstyle 2}
  \right\|^2
  -
  \left\|
  \operator{E}(t)\,\frac{\scriptstyle\phi-\psi}{\scriptstyle 2}
  \right\|^2
  +i\,
  \left\|
  \operator{E}(t)\,\frac{\scriptstyle \phi+i\,\psi}{\scriptstyle 2}
  \right\|^2
  -i\,
  \left\|
  \operator{E}(t)\,\frac{\scriptstyle\phi-i\,\psi}{\scriptstyle 2}
  \right\|^2
\end{equation*}
to define Darboux--Stieltjes integral with respect to
$\rho(t)=\langle \phi,\operator{E}(t)\,\psi\rangle$ as a linear
combination of four Darboux--Stieltjes integrals. Finally, if $f$ is
Darboux--Stieltjes integrable with respect to $\|\operator{E}(t)\,\phi\|^2$
for all $\phi$ in the Hilbert space, then we say that $f$ is Darboux--Stieltjes
integrable with respect to $\operator{E}(t)$.

We aim to extend the convergence of Theorem~\ref{thm:bounded-function-srs}
for continuous functions to Darboux-- or Riemann--Stieltjes
integrable functions.
We start by proving equivalent of Theorem~\ref{thm:bounded-function-srs}
for step functions.
\begin{lemma}
  \label{lemma:step-function-conv}
  Let $\operator{A}_n$ and $\operator{A}$ be self-adjoint operators with
  spectral families $\operator{E}_n(t)$ and $\operator{E}(t)$, respectively, for $n=0,1,2,\ldots$.
  If $\operator{A}_n\xrightarrow{srs}\operator{A}$, then
  for a step function
  $f(\operator{A}_n)\xrightarrow{s}f(\operator{A})$ provided that
  discontinuity points of $f$ are not eigenvalues of $\operator{A}$.
\end{lemma}
\begin{proof}
  A step function can be represented as a finite linear combination of
  characteristic functions. By Theorem~\ref{thm:characteristic-convergence},
  the convergence is strong for each characteristic function. A finite
  linear combination of operators that converge strongly also converges
  strongly \cite[Chapter~4.9, Problem~2]{Kreyszig:1989}.
\end{proof}

Next, we extend the convergence from step functions to
Darboux--Stieltjes integrable functions.
\begin{theorem}
  \label{thm:Darboux-convergence}
  Let $\operator{A}_n$ and $\operator{A}$ be self-adjoint operators with
  spectral families $\operator{E}_n(t)$ and $\operator{E}(t)$, respectively, for $n=0,1,2,\ldots$.
  Let function $f$ be continuous on eigenvalues of $\operator{A}$.
  Let $\operator{A}_n\xrightarrow{srs}\operator{A}$.
  Then the following statements hold:
  \begin{enumerate}
  \item If $f(t)$ is Darboux--Stieltjes integrable with respect to
    $\rho(t)=\|\operator{E}(t)\,\phi\|^2$, then
    \begin{equation*}
      \lim_{n\rightarrow \infty}\langle \phi, f(\operator{A}_n)\,\phi\rangle
      =\langle \phi, f(\operator{A})\,\phi\rangle.
    \end{equation*}
  \item If $f(t)$ is Darboux--Stieltjes integrable with respect to
    $\rho(t)=\langle \psi,\operator{E}(t)\,\phi\rangle$, then
    \begin{equation*}
      \lim_{n\rightarrow \infty}\langle \psi, f(\operator{A}_n)\,\phi\rangle
      =\langle \psi, f(\operator{A})\,\phi\rangle.
    \end{equation*}
  \item If $f(t)$ is Darboux--Stieltjes integrable with respect to
    $\operator{E}(t)$, then
    %
    $
      f(\operator{A}_n)\xrightarrow{s}f(\operator{A}).
$
      %
  \end{enumerate}
\end{theorem}
\begin{proof}
  %
  For the proof of statement 1,
    we take the idea from \cite[Chapter~2.7.8]{Davis+Rabinowitz:1984}.
    We assume that $f$ is Darboux--Stieltjes integrable with respect to
    $\rho(t)=\|\operator{E}(t)\,\phi\|^2$.
    We define step functions $f_l$ and $f_u$ as in
    Lemma~\ref{lemma:Darboux--Stieltjes}. By Lemma~\ref{lemma:Darboux--Stieltjes},
    we can also select $f_l$ and $f_u$ so that they are continuous at
    eigenvalues of $\operator{A}$. By Lemma~\ref{lemma:step-function-conv},
    we have
    \begin{align*}
      f_l(\operator{A}_n)&\xrightarrow{s}f_l(\operator{A}),\\
      f_u(\operator{A}_n)&\xrightarrow{s}f_u(\operator{A}).
    \end{align*}
    Because the strong convergence implies weak convergence, we have also
    \begin{align*}
      \langle \phi,f_l(\operator{A}_n)\,\phi\rangle
      &\rightarrow
      \langle \phi, f_l(\operator{A})\,\phi\rangle,\\
      \langle \phi,f_u(\operator{A}_n)\,\phi\rangle
      &\rightarrow
      \langle \phi, f_u(\operator{A})\,\phi\rangle.
    \end{align*}
    We see from the definition of $f_l$ and $f_u$ in
    Lemma~\ref{lemma:Darboux--Stieltjes} that
    because, for example,
    \begin{align*}
      \langle \phi, f_l(\operator{A})\,\phi\rangle
      &=
      \int_{\sigma(\operator{A})}f_l(t)\,d\|\operator{E}(t)\,\phi\|^2
      \leq
      \int_{\sigma(\operator{A})}f(t)\,d\|\operator{E}(t)\,\phi\|^2
      =
      \langle \phi, f(\operator{A})\,\phi\rangle,
    \end{align*}
    we have
    \begin{equation*}
      \langle \phi, f_l(\operator{A})\,\phi\rangle\leq
      \langle \phi, f(\operator{A})\,\phi\rangle\leq
      \langle \phi, f_u(\operator{A})\,\phi\rangle,
    \end{equation*}
    and
    \begin{equation*}
      \langle \phi, f_l(\operator{A}_n)\,\phi\rangle\leq
      \langle \phi, f(\operator{A}_n)\,\phi\rangle\leq
      \langle \phi, f_u(\operator{A}_n)\,\phi\rangle
    \end{equation*}
    as well.
    By Lemma~\ref{lemma:Darboux--Stieltjes} we also have
    \begin{equation*}
      \langle \phi, (f_u(\operator{A})-f_l(\operator{A}))\,\phi\rangle<\epsilon,
    \end{equation*}
    and thus we get the following inequalities:
    \begin{align*}
      \langle \phi,f(\operator{A})\,\phi\rangle-\epsilon
      &<
      \langle\phi, f_l(\operator{A})\,\phi\rangle
      =
      \liminf_{n\rightarrow\infty} \langle\phi, f_l(\operator{A}_n)\,\phi\rangle,\\
      \langle \phi,f(\operator{A})\,\phi\rangle+\epsilon
      &>
      \langle\phi, f_u(\operator{A})\,\phi\rangle
      =
      \limsup_{n\rightarrow\infty} \langle\phi, f_u(\operator{A}_n)\,\phi\rangle.
    \end{align*}
    By \cite[Theorem~3.17 (b)]{Rudin:1976}, there is N so that for all
    $n\geq N$,
    \begin{equation*}
      \langle \phi,f(\operator{A})\,\phi\rangle-\epsilon
      <
      \langle \phi, f_l(\operator{A}_n)\,\phi\rangle
      \leq
      \langle \phi, f(\operator{A}_n)\,\phi\rangle
      \leq
      \langle \phi, f_u(\operator{A}_n)\,\phi\rangle
      <
      \langle\phi, f(\operator{A})\,\phi\rangle+\epsilon,
    \end{equation*}
    and hence
    \begin{equation*}
      \lim_{n\rightarrow\infty}\langle \phi, f(\operator{A}_n)\,\phi\rangle
      =\langle \phi, f(\operator{A})\,\phi\rangle.
    \end{equation*}

    The proof of statement 2 follows from 1 by the polarization identity.
    For the proof of statement 3, we take an arbitrary $\phi$, and set
    $\psi=f(\operator{A})\,\phi$. Then
    \begin{align*}
      \|(f(\operator{A}_n)-f(\operator{A}))\,\phi\|
      & =
      \langle \phi, f(\operator{A}_n)^2\,\phi\rangle-
      \langle \psi, f(\operator{A}_n)\,\phi\rangle-
      \langle f(\operator{A}_n)\,\phi,\psi\rangle+
      \langle \psi, \psi\rangle\\
      & =
      \langle \phi, f(\operator{A})^2\,\phi\rangle-
      \langle \psi, f(\operator{A})\,\phi\rangle-
      \langle f(\operator{A})\,\phi,\psi\rangle+
      \langle \phi, f(\operator{A})^2\,\phi\rangle\\
      &= 0
    \end{align*}
    as $n\rightarrow\infty$ because $f^2$ is Darboux--Stieltjes integrable
    when $f$ is \cite[Theorem~6.13]{Rudin:1976}.
\end{proof}
By this theorem, we can extend the convergence of
\eqref{eq:bounded-convergence} to functions $f$ that are discontinuous
but Darboux--Stieltjes integrable with respect to $\operator{E}(t)$.

It is good to note that the eigenvalues of $\operator{A}$
are the only points of discontinuity of the family
$\operator{E}(t)$, and points of continuity $t_0$ of
$\operator{E}(t)$ are also the points of continuity for
$\rho(t)=\langle \phi,\operator{E}(t)\,\psi\rangle$.
This is, because by definition,
$\operator{E}(t)\xrightarrow{s}\operator{E}(t_0)$ when
$t\rightarrow t_0$. Strong convergence implies weak convergence
\cite{Riesz+Nagy:1955,Akhiezer+Glazman:1993,Kato:1995,Kreyszig:1989,
Reed+Simon:1981,Weidmann:1980,Simon:2015} which means that
$\langle\phi,\operator{E}(t)\,\psi\rangle\rightarrow
\langle\phi,\operator{E}(t_0)\,\psi\rangle$ for all $\phi,\psi$ in
the Hilbert space. It is still possible that for some
$\phi$ and $\psi$,
$\langle \phi, \operator{E}(t)\,\psi\rangle$ is continuous
at some points that are eigenvalues of $\operator{A}$ and hence
where $\operator{E}(t)$ is not continuous.

If a function $f$
is Riemann--Stieltjes integrable with respect to
$\rho(t)=\|\operator{E}(t)\,\phi\|$ for all $\phi$, then $f$ must be continuous
on all eigenvalues of $\operator{A}$ \cite[p.~251]{Kestelman:1960}.
Then we can say that
$f$ is Riemann--Stieltjes integrable with respect to $\operator{E}(t)$
or with respect to $\operator{A}$.
By Theorem~\ref{thm:Darboux-convergence} statement 3,
we always have convergence for Riemann--Stieltjes integrable functions.
Therefore, we notice that Riemann--Stieltjes integrability is a suitable
criterion for the convergence in terms of discontinuity of the function.
In order to have a discontinuity criterion for unbounded functions alike,
we extend the definition of Riemann--Stieltjes integrability in
the following manner.
\begin{definition}
  Function $f$ is Riemann--Stieltjes integrable with respect to
  self-adjoint operator $\operator{A}$ if $f$ is Riemann--Stieltjes
  integrable with respect to spectral family of $\operator{A}$
  over any finite interval with endpoints that are not eigenvalues
  of $\operator{A}$.
\end{definition}
If an unbounded function $f$ is Riemann--Stieltjes integrable with respect
to self-adjoint operator $\operator{A}$,
then for all $\phi\in\vectorspace{D}(\operator{A})$, there exists
$f(\operator{A})\,\phi=\int_{-\infty}^\infty f(t)\,d\operator{E}(t)\,\phi$
as an improper integral that is a limit of Riemann--Stieltjes integrals
over finite interval.

\subsection{Strong resolvent convergence for unbounded functions}
\label{sec:strong-resolvent-convergence-for-unbounded-functions}
We can now show that an unbounded function of an operator converges in strong resolvent sense. However, as discussed below the result of this theorem looks more useful than it actually is.

\begin{theorem}
  \label{thm:unbounded-function-srs}
  Let $\operator{A}_n$ and $\operator{A}$ be self-adjoint operators, and let
  $f$ be such a continuous real function defined on $\mathbb{R}$
  that $f(\operator{A}_n)$ and $f(\operator{A})$ are densely defined.
  Then
  \begin{equation*}
    \operator{A}_n\xrightarrow{srs}\operator{A}
    \Rightarrow
    f(\operator{A}_n)\xrightarrow{srs}f(\operator{A}).
  \end{equation*}
\end{theorem}
\begin{proof}
  Operators $\operator{B}_n=f(\operator{A}_n)$ and
  $\operator{B}=f(\operator{A})$ are self-adjoint.
  Function
  \begin{equation*}
    \tilde{f}(x)=\frac{1}{f(x)-z},
  \end{equation*}
  is bounded and continuous for nonreal $z$.
  By Theorem~\ref{thm:bounded-function-srs},
  $\tilde{f}(\operator{A}_n)\xrightarrow{s}\tilde{f}(\operator{A})$, and
  we have
  %
  $
    (\operator{B}_n-z)^{-1}=\tilde{f}(\operator{A}_n)
    \xrightarrow{s}\tilde{f}(\operator{A})=(\operator{B}-z)^{-1}
$
    %
  which means
  by definition of the strong resolvent convergence that
  $f(\operator{A}_n)=\operator{B}_n\xrightarrow{srs}\operator{B}=
  f(\operator{A})$.
\end{proof}
It is now important to note that although the above theorem ensures
the strong resolvent convergence, it does not guarantee convergence of the inner product
$\langle \phi, f(\operator{A}_n)\,\phi\rangle$ which is what we actually want.
It only guarantees existence of a sequence $\phi_n\rightarrow \phi$ for which
$\langle \phi, f(\operator{A}_n)\,\phi_n\rangle$ converges as is shown
by the following theorem.
\begin{theorem}
  \label{thm:graph-convergence}
  Let self-adjoint operators $\operator{A}_n$ and $\operator{A}$ be
  such that $\operator{A}_n\xrightarrow{srs}\operator{A}$.
  Then
  for each $\phi\in\vectorspace{D}(\operator{A})$ there is a sequence
  $\phi_n\in\vectorspace{D}(\operator{A}_n)$ such that
  \begin{equation*}
    \|\phi_n-\phi\|+\|\operator{A}_n\,\phi_n-\operator{A}\,\phi\|\rightarrow 0
  \end{equation*}
  as $n\rightarrow \infty$.
\end{theorem}
\begin{proof}
  See \cite[Satz (14)]{Stummel:1972}.
\end{proof}

Example sequence $\phi_n$ is given by $\phi_n=(\operator{A}_n-z)^{-1}\,(\operator{A}-z)\,\phi$, where $z$ is any nonreal number. For nonreal $z$, the sequence $\phi_n$ is not real-valued.
However, it is often possible to find a real vector sequence $(\phi_n)$
based on a dominating
function $h^2\geq |f|$ such that
$\langle\phi_n,f(\operator{A}_n)\,\phi_n\rangle\rightarrow\langle\phi,\operator{A}\,\phi\rangle$.
In the case of numerical integration this changes the construction of the
weights of the quadrature rule so that they are based on the function
sequence
$(\phi_n)$ instead of a fixed function $\phi$ but the convergence
of such quadrature can be proved as is shown by the following theorem.
\begin{theorem}
  \label{thm:reweighted-convergence}
  Let self-adjoint operators $\operator{A}_n$ and $\operator{A}$ be
  such that $\operator{A}_n\xrightarrow{srs}\operator{A}$.
  Let $f$ and $h$ be Riemann--Stieltjes integrable functions with respect to
  $\operator{A}$ so that
  $|f(\mat{x})|\leq h^2(\mat{x})$,
  and $h^2(\mat{x})\geq c > 0$ for all $x\in\Omega$.
    For $\phi\in\vectorspace{D}(h(\operator{A}))$
  \begin{equation*}
    (h(\operator{A}_n))^{-1}\,f(\operator{A}_n)\,
    (h(\operator{A}_n))^{-1}\,h(\operator{A})\,\phi
    \rightarrow (h(\operator{A}))^{-1}\,f(\operator{A})\,\phi
  \end{equation*}
  and hence
  \begin{equation*}
    \langle h(\operator{A})\,\phi, (h(\operator{A}_n))^{-1}\,f(\operator{A}_n)\,
    (h(\operator{A}_n))^{-1}\,h(\operator{A})\,\phi\rangle
    \rightarrow \langle \phi,f(\operator{A})\,\phi\rangle.
  \end{equation*}
\end{theorem}
\begin{proof}
  Function $q(\mat{x})=\frac{f(\mat{x})}{h^2(\mat{x})}$ is bounded so that
  $|q(\mat{x})|\leq 1$. Function $q$ is also Riemann--Stieltjes integrable
  \cite[Theorems 6.11 and 6.13]{Rudin:1976}.
  Therefore, operators
  $\operator{B}_n=(h(\operator{A}_n))^{-1}\,f(\operator{A}_n)\,
  (h(\operator{A}_n))^{-1}$
  are bounded,
  operator norm $\|\operator{B}_n\|\leq 1$ and
  \begin{equation*}
    (h(\operator{A}_n))^{-1}\,f(\operator{A}_n)\,
    (h(\operator{A}_n))^{-1}
    \xrightarrow{s}
    (h(\operator{A}))^{-1}\,
    f(\operator{A})\,(h(\operator{A}))^{-1}.
  \end{equation*}
  Because $\phi\in \vectorspace{D}(h(\operator{A}))$, the result follows.
\end{proof}
\begin{remark}
  This theorem suggests us a new form of convergent quadrature rule.
  If we can compute
  the Fourier
  coefficients $v_i=\langle \phi_i, h\rangle$, using the above theorem, we can
  construct a sequence of convergent quadrature rules as
  \begin{equation*}
    \int_\Omega f(g(\mat{x}))\,w(\mat{x})\,d\mat{x}\approx
    \tilde{\mat{v}}_n^*\,f(\mat{M}_n[g])\,\tilde{\mat{v}}_n,
  \end{equation*}
  where vector $\tilde{\mat{v}}_n=(h(\mat{M}_n[g]))^{-1}\,\mat{v}$ and
  components
  of vector $\mat{v}$ are $v_i$ respectively.
  The nodes are then the eigenvalues
  $\lambda_j$ of each finite matrix approximation $\mat{M}_n[g]$ and the weights are
  \begin{equation}
    \label{eq:general-weights}
    w_j=\frac{|\mat{v}^*\,\mat{u}_j|^2}{h^2(\lambda_j)},
  \end{equation}
  where
  $\mat{u}_j$ are eigenvectors of $\mat{M}_n[g]$ corresponding to the eigenvalue
  $\lambda_j$. Numerical integral of \eqref{eq:numint-def} is a special case
  of this approach with $h=1$. We give a numerical example of the effects
  of this theorem
  in practice in \ref{app:numerical-example}.
\end{remark}

\subsection{Convergence for quadratically bounded functions}
\label{sec:convergence-for-quadratically-bounded-functions}
For convergence without an approximating sequence for
vectors $\phi\in\vectorspace{D}(f(\operator{A}))$, we can prove
a limited form of convergence for a limited class of unbounded
functions.

The following theorem shows that if we can show convergence for one function,
convergence will follow for all smaller functions.
\begin{theorem}
  \label{thm:scalar-dominated-convergence}
  Let $\operator{A}_n$ and $\operator{A}$ be self-adjoint operators
  such that $\operator{A}_n\xrightarrow{srs}\operator{A}$. Let
  $f$ and $h$ be such unbounded functions that they are
  Riemann--Stieltjes integrable with respect
  to $\operator{A}$,
  $|f(x)|\leq |h(x)|$ for all $x\in\mathbb{R}$ and
  $f$ and $h$ are continuous at eigenvalues of $\operator{A}$.
  Then
  for any $\phi\in\vectorspace{D}(h(\operator{A}))$
  \begin{equation*}
    \langle \phi, h(\operator{A}_n)\,\phi\rangle\rightarrow
    \langle \phi, h(\operator{A})\,\phi\rangle
    \Rightarrow
    \langle \phi, f(\operator{A}_n)\,\phi\rangle\rightarrow
    \langle \phi, f(\operator{A})\,\phi\rangle.
  \end{equation*}
\end{theorem}
\begin{proof}
  This follows from \cite[Th\'eor\`eme~3]{Jouravsky:1928}.
\end{proof}

The following lemma makes it possible to show not only the convergence of the
norms of the vectors but the vectors themselves as well.
\begin{lemma}
  \label{lemma:convergence-of-norm}
  Let self-adjoint operators $\operator{A}_n$ converge in strong resolvent
  sense to self-adjoint operator $\operator{A}$. Let
  function
  $h:~\mathbb{R}\mapsto\mathbb{C}$ and
  vector
  $\phi\in\vectorspace{D}(h(\operator{A}))
  \cap_{n=0}^\infty\vectorspace{D}(h(\operator{A}_n))$.
  If
  $\|h(\operator{A}_n)\,\phi\|\rightarrow\|h(\operator{A})\,\phi\|$, then
  also $h(\operator{A}_n)\,\phi\rightarrow h(\operator{A})\,\phi$.
\end{lemma}
\begin{proof}
  Let $\operator{E}(t)$ and $\operator{E}_n(t)$ be the spectral families
  of $\operator{A}$ and $\operator{A}_n$ respectively.
  We can select an
  interval $I$ with endpoints that are not eigenvalues of $\operator{A}$
  such that
  \begin{equation}
    \label{eq:norm-in-tail}
    \int_{\mathbb{R}\setminus I} |h(t)|^2\,d\|\operator{E}(t)\,\phi\|^2<
    \frac{\epsilon^2}{16}=\epsilon_1^2.
  \end{equation}
  Because $\operator{A}_n\xrightarrow{srs}\operator{A}$,
  for all $\epsilon > 0$, there is an $N_1$ such that for all $k>N_1$,
  we have
  \begin{equation}
    \label{eq:convergence-of-norm-in-center}
    \left|
    \int_I |h(t)|^2\,d\|\operator{E}(t)\,\phi\|^2-
    \int_I |h(t)|^2\,d\|\operator{E}_k(t)\,\phi\|^2
    \right| < \frac{\epsilon^2}{8}=\epsilon_2^2
  \end{equation}
  and
  \begin{equation}
    \label{eq:convergence-in-center}
    \left\|
    \int_I h(t)\,d\operator{E}(t)\,\phi-\int_Ih(t)\,d\operator{E}_k(t)\,\phi
    \right\|<\frac{\epsilon}{4}=\epsilon_3.
  \end{equation}
  The convergence of
  $\|h(\operator{A}_n)\,\phi\|\rightarrow\|h(\operator{A})\,\phi\|$ means that
  for all $\epsilon > 0$, there is an $N_2$ such that for all $k>N_2$, we have
  \begin{align*}
    \epsilon_4=\frac{\epsilon^2}{16}
    >&
    |\|h(\operator{A})\,\phi\|^2-\|h(\operator{A}_k)\,\phi\|^2|\\
    =&
    \left|
    \int_I|h(t)|^2\,d\|\operator{E}(t)\,\phi\|^2+
    \int_{\mathbb{R}\setminus I}|h(t)|^2\,d\|\operator{E}(t)\,\phi\|^2
    \right.\\
    &
    \left.
    -\int_I|h(t)|^2\,d\|\operator{E}_k(t)\,\phi\|^2-
    \int_{\mathbb{R}\setminus I}|h(t)|^2\,d\|\operator{E}_k(t)\,\phi\|^2
    \right|,
  \end{align*}
which means by \eqref{eq:norm-in-tail} and
\eqref{eq:convergence-of-norm-in-center} that
\begin{equation}
  \label{eq:convergence-of-norm-in-tail}
  \int_{\mathbb{R}\setminus I}|h(t)|^2\,d\|\operator{E}_k(t)\,\phi\|^2<
  \epsilon_1^2+\epsilon_2^2 +\epsilon_4^2 =\frac{\epsilon^2}{4}=\epsilon_5,
\end{equation}
for all $k>M=\max(N_1,N_2)$.
By \eqref{eq:norm-in-tail}, \eqref{eq:convergence-in-center} and
  \eqref{eq:convergence-of-norm-in-tail} we have
\begin{align*}
  &\|(h(\operator{A})-h(\operator{A}_k))\,\phi\|\\
  &=
  \left\|
  \int_{\mathbb{R}\setminus I} h(t)\,d\operator{E}(t)\,\phi+
  \int_I h(t)\,d\operator{E}(t)\,\phi
  -\int_{\mathbb{R}\setminus I}h(t)\,d\operator{E}_k(t)\,\phi
  -\int_I h(t)\,d\operator{E}_k(t)\,\phi
  \right\|\\
  &\leq
  \left\|
  \int_Ih(t)\,d\operator{E}(t)\,\phi-\int_Ih(t)\,d\operator{E}_k(t)\,\phi
  \right\|
  +\sqrt{\int_{\mathbb{R}\setminus I}|h(t)|^2\,d\|\operator{E}(t)\,\phi\|^2}\\
  &~~~+\sqrt{\int_{\mathbb{R}\setminus I}|h(t)|^2\,d\|\operator{E}_k(t)\,\phi\|^2}\\
  &<
  \epsilon_3+\epsilon_1+\epsilon_5=\epsilon
\end{align*}
for all $k>M$.
\end{proof}

\begin{theorem}
  \label{thm:dominated-convergence}
  Let $\operator{A}_n$ and $\operator{A}$ be self-adjoint operators
  such that $\operator{A}_n\xrightarrow{srs}\operator{A}$. Let
  $f$ and $h$ be such unbounded functions that they are
  Riemann--Stieltjes integrable with respect
  to $\operator{A}$,
  $|f(x)|\leq |h(x)|$ for all $x\in\mathbb{R}$ and
  $f$ and $h$ are continuous at eigenvalues of $\operator{A}$.
  Then
  for any $\phi\in\vectorspace{D}(h(\operator{A}))$
  \begin{equation*}
    h(\operator{A}_n)\,\phi\rightarrow h(\operator{A})\,\phi
    \Rightarrow
    f(\operator{A}_n)\,\phi\rightarrow f(\operator{A})\,\phi.
  \end{equation*}
\end{theorem}
\begin{proof}
This follows from Theorem~\ref{thm:scalar-dominated-convergence} and
  Lemma~\ref{lemma:convergence-of-norm} by using
  the identity
  $\|f(\operator{A})\,\phi\|^2= \langle \phi, |f(\operator{A})|^2\,\phi \rangle$
  for functions $f$ and $h$.
\end{proof}

For example, we can apply this theorem to linearly bounded functions.
\begin{theorem}
  \label{thm:linearly-bounded-convergence}
  Let operator $\operator{A}$ be self-adjoint with an infinite matrix
  representation $\mat{A}_{\infty}$ and vector $\mat{v}\in\ell^2_0$.
  Let finite approximations be
  \begin{align*}
    \left[
      \mat{A}_n
      \right]_{i,j}
    &=
    \left[
      \mat{A}_\infty
      \right]_{i,j},\\
    [\mat{v}_n]_i
    &=
    [\mat{v}]_i
  \end{align*}
  for $i,j=0,1,\ldots,n$.
  Let function $f$ be Riemann--Stieltjes integrable with respect to
  $\operator{A}$ and
  linearly bounded, that is, for all $x\in\mathbb{R}$,
  $|f(x)|\leq a + b|x|$
  for some positive $a$ and $b$.
  Then
  \begin{equation*}
    f(\mat{A}_n)\,\mat{v}_n
    \rightarrow
    f
    \left(
    \mat{A}_\infty
    \right)
    \,\mat{v}.
  \end{equation*}
\end{theorem}
\begin{proof}
  By Theorem~\ref{thm:srs-convergence},
  the convergence holds for function $h(x) = x$.
  By Theorem~\ref{thm:dominated-convergence},
  the convergence
  holds also for function $h(x)=|x|$. The convergence then holds
  also for function $h(x) = b\,|x|$, for $h(x) = a+b\,|x|$
  and again by Theorem~\ref{thm:dominated-convergence} for $f(x)$.
\end{proof}
We can apply this theorem to a product of two linearly bounded
functions.
\begin{theorem}
  \label{thm:product-convergence}
  Let self-adjoint operators $\operator{A}$ and $\operator{B}$ have
  infinite matrix representations $\mat{A}_\infty$ and $\mat{B}_\infty$.
  Furthermore, let vectors $\mat{u},\mat{v}\in\ell^2_0$ and
  finite approximations have the elements
  \begin{align*}
    [\mat{A}_n]_{i,j}&=[\mat{A}_\infty]_{i,j},\\
    [\mat{B}_n]_{i,j}&=[\mat{B}_\infty]_{i,j},\\
    [\mat{u}_n]_i  &= [\mat{u}]_i,\\
    [\mat{v}_n]_i &= [\mat{v}]_i,
  \end{align*}
  for $i,j=0,1,\ldots,n$.
  Then for linearly bounded $f_1,f_2$ that are Riemann--Stieltjes
  integrable with respect to $\operator{A}$ and $\operator{B}$ respectively,
  we have
    \begin{equation*}
      \mat{u}_n^*\,
      f_1\left(
      \mat{A}_n
      \right)\,
      f_2\left(
      \mat{B}_n
      \right)\,
      \mat{v}_n
      \rightarrow
      \mat{u}^*\,f_1(\mat{A}_\infty)\,f_2(\mat{B}_\infty)\,\mat{v}.
    \end{equation*}
\end{theorem}
\begin{proof}
  This follows from
  Theorem~\ref{thm:linearly-bounded-convergence} and
  the property that in a Hilbert space, if
  $\mat{x}_n\rightarrow \mat{x}$ and $\mat{y}_n\rightarrow \mat{y}$ then
  $\langle \mat{x}_n,\mat{y}_n\rangle \rightarrow \langle \mat{x},\mat{y}\rangle$
  \cite[p.~199]{Riesz+Nagy:1955} or \cite[Lemma~3.2-2]{Kreyszig:1989}.
\end{proof}

We can now apply this theorem to quadratically bounded functions.
\begin{theorem}
  \label{thm:quadratic-convergence}
  Let self-adjoint operator $\operator{A}$ have
  infinite matrix representation $\mat{A}_\infty$.
  Let vectors $\mat{u},\mat{v}\in\ell^2_0$, and
  let finite approximations have the
  elements
  \begin{align*}
    [\mat{A}_n]_{i,j}&=[\mat{A}_\infty]_{i,j},\\
    [\mat{u}_n]_i  &= [\mat{u}]_i,\\
    [\mat{v}_n]_i &= [\mat{v}]_i,
  \end{align*}
  for $i,j=0,1,\ldots,n$.
  Let function $f$ be Riemann--Stieltjes integrable
  with respect to $\operator{A}$
  and quadratically bounded, that is,
  for all $x\in\mathbb{R}$,
    $|f(x)|\leq a + b\,x^2$, for some positive $a$ and $b$. Then
    \begin{equation*}
      \mat{u}_n^*\,
      f\left(
      \mat{A}_n
      \right)
      \,
      \mat{v}_n
      \rightarrow
      \mat{u}^*\,f(\mat{A}_\infty)\,\mat{v}.
    \end{equation*}
\end{theorem}
\begin{proof}
  This follows from Theorem~\ref{thm:product-convergence} by selecting
  $\operator{A}=\operator{B}$,
  $f_1=\operatorname{sgn}(f)\,\sqrt{|f|}$ and $f_2=\sqrt{|f|}$, where
  $\operatorname{sgn}(x)$
  is the sign function. The product function $f_1\,f_2$ is Riemann--Stieltjes
  integrable by \cite[Theorem 6.11 and 6.13]{Rudin:1976}
\end{proof}

For the convergence of the numerical approximation of integrals, we
then have the following theorem.
\begin{theorem}
  Let orthonormal functions $\phi_0=1,\phi_1,\phi_2,\ldots$ be dense in
  $\vectorspace{L}^2_w(\Omega)$.
  Let function
  $g$ satisfy \eqref{eq:mult-mat} and
  let matrices $\mat{M}_n[g]$ have
  elements as in \eqref{eq:matrix-elements}.
  Further, let quadratically bounded function $f$ be
  Riemann--Stieltjes integrable with respect to $\operator{M}[g]$.
  Then
  \begin{equation*}
    \lim_{n\rightarrow\infty} [f(\mat{M}_n[g])]_{0,0}=
    \int_{\Omega} f(g(\mat{x}))\,w(\mat{x})\,d\mat{x}.
  \end{equation*}
  Similarly, for functions $g_1$ and $g_2$ satisfying \eqref{eq:mult-mat}
  and linearly bounded $f_1$ and $f_2$ we have
  \begin{equation*}
    \lim_{n\rightarrow\infty} [f_1(\mat{M}_n[g_1])\,f_2(\mat{M}_n[g_2])]_{0,0}=
    \int_{\Omega} f_1(g_1(\mat{x}))\,f_2(g_2(\mat{x}))\,w(\mat{x})\,d\mat{x}.
  \end{equation*}
\end{theorem}
\begin{proof}
The result follows from Theorems
\ref{thm:product-convergence}, and \ref{thm:quadratic-convergence}.
\end{proof}
We can extend this result to polynomially bounded functions in the
special case that the inside function $g$ and the basis
functions are polynomials.
We define the total degree of a multivariate monomial
$\prod_{i=0}^d x_i^{n_i}$ as $\sum_{i=0}^d n_i$. Then for a polynomial,
that is, a finite linear combination of monomials, its total degree is the
total degree of the highest monomial.
If polynomials are dense in a Hilbert space, they can
be partially ordered according to the total degree of the polynomials,
for example, in graded lexicographic order
\cite[Chapter~3.1]{Dunkl+Xu:2014}.
\begin{theorem}
  Let the orthonormal polynomials $\phi_0,\phi_1,\phi_2,\ldots$ be partially
  ordered by their total degree and dense in
  $\vectorspace{L}^2_w(\Omega)$.
  Let the function $g$ be polynomial and
  let the matrices $\mat{M}_n[g]$ have
  elements as in \eqref{eq:matrix-elements}. Let a
  function $f$ be polynomially bounded
  that is,
  $f(x)\leq a + b\,|x|^m$ for some positive $a,b$ and $m\in\mathbb{N}$, then
  \begin{equation*}
    \lim_{n\rightarrow\infty} [f(\mat{M}_n[g])]_{0,0}=
    \int_{\Omega} f(g(\mat{x}))\,w(\mat{x})\,d\mat{x},
  \end{equation*}
  if $f$ is  Riemann--Stieltjes integrable with respect to $\operator{M}[g]$.
\end{theorem}
\begin{proof}
  Because polynomials are dense in   $\vectorspace{L}^2_w(\Omega)$, the function
  $a+b\,|g(\mat{x})|^m$ is integrable and $g$ satisfies \eqref{eq:mult-mat}.
  Let the total degree of polynomial $g$ be $k$.
  Let the projection operator $\operator{P}_n$ be defined by
  \begin{equation*}
    \operator{P}_n\,\psi=\sum_{i=0}^n\langle \phi_i,\psi\rangle\,\phi_i.
  \end{equation*}
  For a polynomial $\phi_j$ of total degree $p$, we can
  select $n$ so that the linear combination of
  $\{\phi_i\}_{i=0}^n$ covers all polynomials of
  total degree $k\,m + p$ and then we have the following equalities
  \begin{align*}
    g^m\,\phi_j = \operator{M}[g]^m\,\phi_j
    &=
    \left(
    \operator{P}_n\,\operator{M}[g]\,\operator{P}_n
    \right)^m\,\phi_j\\
    \mat{M}_\infty[g]^m\,\mat{e}_j
    &=
    \left[
      \begin{array}{ccc}
        \mat{M}_n[g] & 0      & \ldots \\
        0            & 0      & \ldots \\
        \vdots       & \vdots & \ddots
      \end{array}
      \right]^m\,\mat{e}_j.
  \end{align*}
  By Theorem~\ref{thm:dominated-convergence} we have
  \begin{equation*}
    \left[
      \begin{array}{ccc}
        f(\mat{M}_n[g]) & 0      & \ldots \\
        0            & 0      & \ldots \\
        \vdots       & \vdots & \ddots
      \end{array}
      \right]\,\mat{e}_j
    \rightarrow
      f(\mat{M}_\infty[g])\,\mat{e}_j,
  \end{equation*}
  from where it follows that
  \begin{equation*}
  [f(\mat{M}_n[g])]_{0,0}=\mat{e}_0^\top\,f(\mat{M}_n[g])\,\mat{e}_0\rightarrow
  \mat{e}_0^\top\,f(\mat{M}_\infty[g])\,\mat{e}_o=
  \int_\Omega f(g(\mat{x}))\,w(\mat{x})\,d\mat{x}.
  \end{equation*}
\end{proof}
\begin{remark}
  The same arguments hold more generally for convergence of
  $\mat{e}_i^\top\,f(\mat{A}_n)\,\mat{e}_j\rightarrow \mat{e}_i^\top\,f(\mat{A}_\infty)\,\mat{e}_j$
  when $\mat{A}_\infty$ is a self-adjoint infinite sparse matrix, that is,
  an infinite matrix with only finite number of non-zero components in each
  column and row.
\end{remark}
\begin{remark}
  In the special case when the inside function is
  $g(x)=\operatorname{id}(x)=x$, we have essentially the same proof as
  in \cite[Sections~1.8~and~1.13]{Simon:2009} for the convergence of
  Gaussian quadrature on
  polynomially bounded functions.
\end{remark}

We can use the following lemma to prove convergence beyond polynomially bounded
functions.
\begin{lemma}
  \label{lemma:bounded-convergence}
  Let self-adjoint operators $\operator{A}_n$ converge to a self-adjoint
  operator $\operator{A}$ in strong resolvent sense. Let
  function
  $f:~\mathbb{R}\mapsto [0,\infty)$ and vector $\phi$ be such that
  $\langle \phi,f(\operator{A})\,\phi\rangle<\infty$.
  If
  $\langle \phi, f(\operator{A}_n)\,\phi\rangle$ is bounded by
  $\langle \phi,f(\operator{A})\,\phi\rangle$, then
  $\langle \phi, f(\operator{A}_n)\,\phi\rangle \rightarrow
  \langle \phi, f(\operator{A})\,\phi\rangle$.
\end{lemma}
\begin{proof}
Let the spectral families of $\operator{A}_n$ and $\operator{A}$ be
$\operator{E}_n(t)$ and $\operator{E}(t)$ respectively.
Because $f$ is nonnegative, for an $\epsilon > 0$, we can select a finite
interval $I$ with endpoints that are not eigenvalues of $\operator{A}$, so that
\begin{equation*}
  \langle \phi, f(\operator{A})\,\phi\rangle - \frac{\epsilon}{2} <
  \int_I f(t)\,d\langle \phi,\operator{E}(t)\,\phi\rangle \leq
  \langle \phi, f(\operator{A})\,\phi\rangle.
\end{equation*}
Due to strong resolvent convergence, there is
$N\in\mathbb{N}$, so that for all $k>N$, we have
\begin{equation*}
  \left|
  \int_I f(t)\,d\langle \phi, \operator{E}(t)\,\phi\rangle -
  \int_I f(t)\,d\langle \phi, \operator{E}_k(t)\,\phi \rangle
  \right| < \frac{\epsilon}{2}.
\end{equation*}
Because, $f$ is nonnegative, these inequalities mean that for all $k>N$,
we have
\begin{equation*}
  \langle \phi, f(\operator{A})\,\phi\rangle - \epsilon
  <
  \int_I f(t)\,d\langle \phi, \operator{E}_k(t)\,\phi\rangle
  \leq
  \langle \phi, f(\operator{A}_k)\,\phi \rangle
  \leq
  \langle \phi, f(\operator{A})\,\phi \rangle.
\end{equation*}
\end{proof}

In \cite[Chapter~IV, Section~8--10]{Shohat+Tamarkin:1963},
the proof of convergence for
  Gaussian quadratures is based on the inequality
  \begin{equation}
    \label{eq:gauss-inequality}
    \sum_{j=0}^n w_j\,h(x_j) \leq \int_\Omega h(x)\,w(x)\,dx,
  \end{equation}
  when all even derivatives of $h$ are nonnegative at origin which follows from
  \eqref{eq:wx_ineq}. This kind of
  inequality makes it possible to prove convergence for functions that are
  bounded by functions that can have negative odd derivatives.
  When the basis functions are polynomials,
  in terms of the infinite matrices and the leading principal submatrices,
  \eqref{eq:gauss-inequality} can be expressed as
  \begin{equation}
    \label{eq:matrix-gauss-inequality}
    [h(\mat{M}_n[\operatorname{id}])]_{0,0} \leq
    [h(\mat{M}_\infty[\operatorname{id}])]_{0,0},
  \end{equation}
  where function $\operatorname{id}(x)=x$ but it does not hold for
  arbitrary inside functions $g\neq \operatorname{id}$ when the basis
  functions are
  not polynomials or for
 other matrix elements except
  the $0,0$ element. Therefore, we can consider
  functions with nonnegative even derivatives at zero the highest possible
  bound for which we have convergence
  if we can generalize \eqref{eq:matrix-gauss-inequality} for other
  elements than just $0,0$ element and arbitrary
  basis functions and inside function $g$.

\subsection{Nonnegative matrix coefficients}
\label{sec:nonnegative-matrix-coefficients}
We can further expand the space of numerically integrable functions
when the matrix coefficients are all nonnegative. This can seem like a hard
restriction but, for example, many Jacobi matrices are such.
For example, \cite[Table~1.1]{Gautschi:2004} lists Jacobi matrix coefficients
for about dozen different so called classical orthogonal polynomials.
The only ones with possibly negative Jacobi matrix coefficients are the
Chebyshev polynomials of fourth kind and the Jacobi and the Meixner--Pollaczek
polynomials
with parameter values on certain ranges.

In some cases, it is also possible to extend the results for matrices with nonnegative coefficients
to matrices with positive and negative real coefficients by considering the positive
and negative parts separately.

We extend the concept of $m$-positive and $m$-monotone functions, so that
matrix entries are allowed to be positive and $a\leq\infty$ for $a<\infty$
but $\infty \leq \infty$ is not true.
\begin{definition}
  \label{def:unbounded-m-positive}
  Let $f:I\mapsto \mathbb{R}$ be a real function defined on an interval $I\subset \mathbb{R}$.
  \begin{enumerate}
  \item
    $f$ is \textit{unbounded $m$-positive} if $0\in I$ and $f(\mat{A}){\scriptstyle\overset{ew}\geq} 0$ for every
    symmetric infinite matrix $\mat{A}{\scriptstyle\overset{ew}\geq} \mat{0}$ with spectrum in $I$.
    $\mat{A}$ and $f(\mat{A})$ are not required to satisfy
    \eqref{eq:infinite-matrix-condition}.
    Elements of
    $f(\mat{A})$ are allowed to be $\infty$.
  \item
    $f$ is \textit{unbounded $m$-monotone} if
    $\mat{A}{\scriptstyle\overset{ew}\leq} \mat{B} \Rightarrow f(\mat{A}){\scriptstyle\overset{ew}\leq} f(\mat{B})$, for all
    self-adjoint infinite matrices $\mat{A},\mat{B}{\scriptstyle\overset{ew}\geq} \mat{0}$ with spectra in $I$.
    $\mat{A}$, $\mat{B}$, $f(\mat{A})$ and $f(\mat{B})$ are not required
    to satisfy \eqref{eq:infinite-matrix-condition}.
    Elements of
     $f(\mat{B})$ are allowed to be $\infty$.
  \end{enumerate}
\end{definition}
For example, monomials and nonnegative linear combinations of monomials
are unbounded $m$-positive functions.
We can extend Theorem~\ref{thm:m-positive-equal-m-monotone} to unbounded
$m$-positive and unbounded $m$-monotone functions.
\begin{theorem}
  \label{thm:operator-m}
  Let $f:I\mapsto \mathbb{R}$ be a real function defined on an interval
  $I\subset \mathbb{R}$ such that $0\in I$ and $f(\mat{0}){\scriptstyle\overset{ew}\geq} \mat{0}$.
  Then $f$ is
  unbounded $m$-positive if and only if it is unbounded $m$-monotone.
\end{theorem}
\begin{proof}
  We adapt the proof of \cite[Theorem~2.1]{Hansen:1992} for symmetric arbitrary
  sized $n\times n$ matrices to unbounded symmetric infinite matrices.
  First, we prove that if $f$ is unbounded $m$-monotone, then it is also
  unbounded $m$-positive. We set $\mat{A}=\mat{0}$ to obtain
  $\mat{0}{\scriptstyle\overset{ew}\leq} f(\mat{0}){\scriptstyle\overset{ew}\leq} f(\mat{B})$ for all symmetric
  infinite matrices.

  To prove that unbounded $m$-positive $f$ is also unbounded $m$-monotone,
  we consider symmetric infinite matrices $\mat{A},\mat{B}{\scriptstyle\overset{ew}\geq} \mat{0}$
  with the spectra in $I$.
  We construct a $2\times 2$ block infinite matrix
  \begin{equation*}
    X=
    \left[
      \begin{array}{cc}
        \mat{A} & \mat{0} \\
        \mat{0} & \mat{B}
      \end{array}
      \right].
  \end{equation*}
  We can extend the ${\scriptstyle\overset{ew}\leq}$ relation for the block matrix simply as blockwise.
  The spectrum of $\mat{X}$ is also contained in $I$.
  We define a $2\times 2$ unitary infinite block matrix
  \begin{equation*}
    U=
    \frac{1}{\sqrt{2}}
    \left[
      \begin{array}{cc}
        I_\infty & -I_\infty\\
        I_\infty & I_\infty
      \end{array}
      \right].
  \end{equation*}
  We calculate
  \begin{equation*}
    \mat{U}^*\,\mat{X}\,\mat{U}=
    \frac{1}{2}
    \left[
      \begin{array}{cc}
        \mat{A}+\mat{B} & \mat{B} - \mat{A}\\
        \mat{B} - \mat{A} & \mat{A} + \mat{B}
      \end{array}
      \right]
  \end{equation*}
  and notice that the spectrum of
  $\mat{U}^*\,\mat{X}\,\mat{U}$ is
  contained in $I$ and that
  $\mat{U}^*\,\mat{X}\,\mat{U}{\scriptstyle\overset{ew}\geq} \mat{0}$
  if $\mat{A}{\scriptstyle\overset{ew}\leq}\mat{B}$. If $f$ is unbounded $m$-positive, then
  \begin{align*}
    \mat{0} &{\scriptstyle\overset{ew}\leq} f(\mat{U}^*\,\mat{X}\,\mat{U})
    =
    \mat{U}^*\,f(\mat{X})\,\mat{U}\\
    &=
    \mat{U}^*\,
    \left[
      \begin{array}{cc}
        f(\mat{A}) & 0\\
        0          & f(\mat{B})
      \end{array}
      \right]
    \mat{U}\\
    &=
    \frac{1}{2}
    \left[
      \begin{array}{cc}
        f(\mat{A})+f(\mat{B}) & f(\mat{B})-f(\mat{A})\\
        f(\mat{B})-f(\mat{A}) & f(\mat{A})+f(\mat{B})
      \end{array}
      \right],
  \end{align*}
  which implies that $f(\mat{A}) {\scriptstyle\overset{ew}\leq} f(\mat{B})$. Here we multiply only
  $f(\mat{A})$, which has finite elements, with $-\mat{I}_\infty$ and define
  $\infty-a=\infty$ for all $a\in[0,\infty)$.
\end{proof}

 %
 %
  %
   %
    %
 %
  %
  %
  %
 %
  %
   %
    %
 %
  %
   %
  %
  %

Absolutely monotone functions, that is, functions, that have power series
expansion $h(x)=\sum_{n=0}^\infty c_n\,x^n$ where $c_n\geq 0$,
are unbounded $m$-positive and
therefore also unbounded $m$-monotone. These functions are also very fast
growing, so if we have convergence for such function, then we have also
convergence for large family of functions that grow slower.
  Examples of functions that have power series representation with
  nonnegative coefficients
  are $\exp(x)$, $\cosh(\sqrt{x})$ that were used in \cite{Uspensky:1928}
  and more generally the Mittag--Leffler functions
  \begin{equation*}
    E_\gamma(x)=\sum_{k=0}^\infty \frac{x^k}{\Gamma(\gamma\,k+1)},
  \end{equation*}
  where $0<\gamma\leq 2$ and $x\in[0,\infty)$
  \cite[eq. (3.3)]{Bultheel+Diaz-Mendoza+Gonzalez-Vera+Orive:2000}.

For functions bounded by unbounded $m$-positive functions we have the following
theorem.
\begin{theorem}
  \label{thm:m-positive-convergence}
  Let a self-adjoint infinite matrix $\mat{A}_\infty{\scriptstyle\overset{ew}\geq} \mat{0}$ satisfy
  \eqref{eq:infinite-matrix-condition}.
  Let function $f$ be bounded by a positive and unbounded $m$-positive
  function $h$, that is, $|f(x)|\leq h(x)$ for $x\in\sigma(\mat{A}_\infty)$.
  Let finite matrix approximations $\mat{A}_n$ be leading principal
  submatrices of $\mat{A}_\infty$.
  Then
  $[f(\mat{A}_n)]_{k,k}\rightarrow [f(\mat{A}_\infty)]_{k,k}$
  for all $k=0,1,2,\ldots$ for which $[h(\mat{A})]_{k,k}<\infty$.
\end{theorem}
\begin{proof}
  Let self-adjoint infinite matrices
  \begin{equation*}
    \mat{B}_n=
    \left[
      \begin{array}{ccc}
        \mat{A}_n & 0      &\ldots\\
        0         & 0      & \ldots\\
        \vdots    & \vdots & \ddots
      \end{array}
      \right].
  \end{equation*}
  Because $\mat{B}_n{\scriptstyle\overset{ew}\leq} \mat{A}$,
  by Theorem~\ref{thm:operator-m}, the sequence
  $[h(\mat{B}_n)]_{k,k}$ is bounded by $[h(\mat{A})]_{k,k}$ which means by
  Lemma~\ref{lemma:bounded-convergence} that it must converge to
  $[h(\mat{A})]_{k,k}$.
  Because for $n\geq k$, we have
  $[h(\mat{A}_n)]_{k,k}=[h(\mat{B}_n)]_{k,k}$, we also have
  $[h(\mat{A}_n)]_{k,k}\rightarrow [h(\mat{A})]_{k,k}$.
  Due to Theorem~\ref{thm:scalar-dominated-convergence}, we have
  convergence for $f$ bounded by $h$.
\end{proof}

In terms of restrictions for the bounds of the outside function $h$,
we can reach the
most generic level of convergence, that is, convergence for functions
that have nonnegative even derivatives, like in
\cite[Chapter~IV, Section~8--10]{Shohat+Tamarkin:1963},
on interval $(-\infty,\infty)$
when $w(x)$ is an even function, basis functions are even or odd
and $g(x)$ is odd. In more general cases,
we have to use a stricter condition that the bounding function $h$
is absolutely monotone, that is, has all derivatives positive.
\begin{theorem}
  \label{thm:power-series-convergence}
  Let the orthonormal functions $\phi_0,\phi_1,\phi_2,\ldots$ be dense in
  $\vectorspace{L}^2_w(\Omega)$.
  Let the function
  $g$ satisfy \eqref{eq:mult-mat} and let the matrices $\mat{M}_n[g]$ have
  elements as in \eqref{eq:matrix-elements} that are nonnegative.
  Furthermore, let a function $f$ be Riemann--Stieltjes integrable with
  respect to $\operator{M}[g]$ and bounded by
  a function $h$ so that $|f(x)|\leq h(x)$ and
    \begin{equation*}
      \int_{\Omega} h(g(\mat{x}))\,|\phi_i(x)|^2\,w(\mat{x})\,d\mat{x} < \infty.
    \end{equation*}
  Then
  \begin{equation*}
    \lim_{n\rightarrow\infty} [f(\mat{M}_n[g])]_{i,i}=
    \int_{\Omega} f(g(\mat{x}))\,|\phi_i(x)|^2\,w(\mat{x})\,d\mat{x},
  \end{equation*}
  if one of the following additional conditions holds:
  \begin{enumerate}
  \item Function $h(x)$ has the form
    \begin{equation*}
      h(x) = \sum_{k=0}^\infty a_k\,x^k,
    \end{equation*}
    where $a_k\geq 0$ for all $k=0,1,\ldots$.
  \item Function $g(\mat{x})$ satisfies
    \begin{equation*}
      \int_\Omega g(\mat{x})^{2\,k+1}\,|\phi_i(x)|^2\,w(\mat{x})\,d\mat{x}
      =[\mat{M}_n[g]^{2\,k+1}]_{i,i}=0
    \end{equation*}
    for all $k=0,1,2,\ldots$ and all even coefficients of
    \begin{equation*}
      h(x)=\sum_{k=0}^\infty a_k\,x^k
    \end{equation*}
    %
    are nonnegative: $a_{2k}\geq 0$.
  \end{enumerate}
\end{theorem}
\begin{proof}
  Proof follows from
  Theorem~\ref{thm:m-positive-convergence}. Under the second condition
  the odd $a_{2k+1}$ coefficients can also be negative.
\end{proof}

We give a numerical example demonstrating the difference in convergence
for matrices with purely nonnegative elements and matrices with
mixed sign elements in \ref{app:numerical-example}.

In some cases, it is also possible to prove convergence when
matrix coefficients have mixed signs by considering the negative and positive
parts separately. We decompose a self-adjoint real infinite matrix $\mat{A}$
into positive $\mat{P}$
and negative $\mat{N}$ parts so that
\begin{align*}
  [\mat{P}]_{i,j}&=
  \left\{
    \begin{array}{llll}
      [\mat{A}]_{i,j} &\mathrm{for} & [\mat{A}]_{i,j}\geq 0,\\
      0              &\mathrm{for} & [\mat{A}]_{i,j}<0,
    \end{array}
    \right.\\
    [\mat{N}]_{i,j}&=
    \left\{
      \begin{array}{lll}
      -[\mat{A}]_{i,j} &\mathrm{for} & [\mat{A}]_{i,j}\leq 0,\\
      0             &\mathrm{for} & [\mat{A}]_{i,j}>0.
      \end{array}
      \right.
\end{align*}
Then $\mat{A}=\mat{P} - \mat{N}$. We construct
a $2\times 2$ block infinite matrix with nonnegative elements
\begin{equation*}
  \left[
    \begin{array}{cc}
      \mat{P} & \mat{N}\\
      \mat{N} & \mat{P}
    \end{array}
    \right].
\end{equation*}
We can express any
matrix element of a function of a self-adjoint real infinite matrix as a
function of a $2\times 2$ block infinite matrix with nonnegative coefficients
as
\begin{equation*}
  [f(\mat{A})]_{i,j}=
  \mat{e}_i^\top\,f(\mat{A})\,\mat{e}_j=
  \left[
    \begin{array}{cc}
      \mat{e}_i & \mat{0}
      \end{array}
    \right]^\top\,
  f\left(\left[
    \begin{array}{cc}
      \mat{P} & \mat{N}\\
      \mat{N} & \mat{P}
    \end{array}
    \right]
  \right)
  \,\left[
    \begin{array}{c}
      \mat{e}_j\\
      -\mat{e}_j
    \end{array}
    \right].
\end{equation*}

For example,
we consider an unbounded function $g(x)=x\,\cos(\sqrt{x})$, a weight
function $w(x)=e^{-x}$ on interval $[0,\infty)$ and polynomials as orthogonal
functions.
  The infinite matrix $\mat{M}_\infty[g]$ is unbounded because function
  $g$ is an unbounded function on interval $[0,\infty)$. The elements of
    $\mat{M}_\infty[g]$ are given by
\begin{align*}
  \mat{M}_\infty[g]
  =&
  \frac{\mathrm{i}\sqrt{\pi}\,\mathrm{erf}
    \left(
    \frac{\mathrm{i}}
         {2}
         \right)
         \,e^{-\frac{1}{4}}}
       {8192}\,
  \left[
    \begin{array}{cccc}
      5120   & 5376   & -1120  & \ldots \\
      5376   & 13632  & 4456   & \ldots \\
      -1120  & 4456   & 10293  & \ldots \\
      \vdots & \vdots & \vdots & \ddots
    \end{array}
    \right]\\
    & +
  \frac{1}{4096}\,
  \left[
    \begin{array}{cccc}
      3072   & 768    & -2208  & \ldots\\
      768    & 192    & -6312  & \ldots\\
      -2208  & -6312  & -11613 & \ldots \\
      \vdots & \vdots & \vdots & \ddots
    \end{array}
    \right],
\end{align*}
where $\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}dt$.
Some of the elements are negative and some positive.
The pattern for the signs of the matrix elements is
\begin{equation*}
  \left[
    \begin{array}{ccccccc}
      + & - & - & + & - & + & \ldots\\
      - & - & - & - & + & - & \ldots\\
      - & - & - & - & + & + & \ldots\\
      + & - & - & - & - & + & \ldots\\
      - & + & + & - & - & - & \ldots\\
      + & - & + & + & - & - & \ldots\\
      \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots
    \end{array}
    \right].
\end{equation*}
However, it is possible to prove the convergence for polynomially bounded
functions as is shown
in the following theorem.
\begin{theorem}
  Let $w(x)=e^{-x}$ on interval $[0,\infty)$,
    the orthogonal basis functions $\phi_i$ polynomials
and $g(x)=x\,\cos(\sqrt{x})$. Let $\mat{M}_n[g]$ have
elements as in  \eqref{eq:matrix-elements}. Then
\begin{equation*}
  [f(\mat{M}_n[g])]_{i,j}\rightarrow [\mat{M}_\infty[f(g)]]_{i,j}=
  \int_0^\infty \phi_i(x)\,\phi_j(x)\,f(g(x))\,w(x)\,dx,
\end{equation*}
when $f$ is bounded by some polynomial.
\end{theorem}
\begin{proof}
The orthogonal polynomials $\phi_i$ are known as the Laguerre polynomials.
Elements of the Jacobi matrix, the infinite matrix representation for
the function $\operatorname{id}(x)=x$, are
nonnegative. We can see this by looking at the
exact formula for the recurrence coefficients of the Jacobi matrix
for the Laguerre polynomials
  \cite[Table~1.1]{Gautschi:2004} or by looking at the few first matrix
  coefficients
  \begin{equation}
    \label{eq:Laguerre-system-Jacobi-matrix}
    \mat{M}_4[\operatorname{id}]=
    \left[
      \begin{array}{ccccc}
        1  &   1  &   0  &   0  &   0\\
        1  &   3  &   2  &   0  &   0\\
        0  &   2  &   5  &   3  &   0\\
        0  &   0  &   3  &   7  &   4\\
        0  &   0  &   0  &   4  &   9
      \end{array}
      \right],
  \end{equation}
  where we see that the main diagonal elements are the odd numbers and
  the off-diagonal
  elements are natural numbers all in a growing order and thus nonnegative.

  For a function $h(x)=x\,\cosh(\sqrt{x})$, the elements of an infinite
  matrix $\mat{M}_\infty[h]$ are also nonnegative because $h$
  has a Taylor series expansion
  \begin{equation*}
    h(x)=x\,\cosh(\sqrt{x})=\sum_{k=0}^\infty \frac{1}{(2\,k)!}\,x^{k+1}
  \end{equation*}
  with only nonnegative coefficients
  and that converges for all $x\in[0,\infty)$.
    The infinite matrix representation of the multiplication operator
    of function $h$ has series representation for the elements as
    \begin{align*}
      [\mat{M}[h]]_{i,j}
      &=
      \sum_{k=0}^\infty \frac{1}{(2\,k)!}\,
          [\mat{M}_\infty[\operatorname{id}]^{k+1}]_{i,j} \geq 0.
    \end{align*}
    The function $g$ has almost the same Taylor series expansion as
    the function $h$
    except some coefficients
    are negative
\begin{equation*}
  g(x)=x\,\cos(\sqrt{x})=\sum_{k=0}^\infty \frac{(-1)^k}{(2\,k)!}\,x^{k+1}.
\end{equation*}
However, the negative coefficients are bounded by the coefficients of
the Taylor series of the function $h$.

We decompose $\mat{M}_\infty[g]$ into positive and negative part as
$\mat{M}_\infty[g]=\mat{P}-\mat{N}$.
The elements of both the positive $\mat{P}$ and
the negative $\mat{N}$ part are bounded by the elements of the infinite matrix
representation
of the multiplication operator for the function $h$ as
\begin{align*}
  [\mat{P}]_{i,j}
  &=
  \sum_{k=0}^\infty \frac{1}{(4\,k)!}\,
      [\mat{M}_\infty[\operatorname{id}]^{2\,k+1}]_{i,j} \leq
      [\mat{M}_\infty[h]]_{i,j},\\
  [\mat{N}]_{i,j}
  &=
  \sum_{k=0}^\infty \frac{1}{(4\,k + 2)!}\,
      [\mat{M}_\infty[\operatorname{id}]^{2\,k+2}]_{i,j} \leq
      [\mat{M}_\infty[h]]_{i,j},
\end{align*}
that is, $\mat{P},\mat{N}{\scriptstyle\overset{ew}\leq} \mat{M}_\infty[h]$.
By Theorem~\ref{thm:operator-m}
\begin{equation*}
  \left[
    \begin{array}{cc}
      \mat{P} & \mat{N}\\
      \mat{N} & \mat{P}
    \end{array}
    \right]^m
  {\scriptstyle\overset{ew}\leq}
  \left[
    \begin{array}{cc}
      \mat{M}_\infty[h] & \mat{M}_\infty[h] \\
      \mat{M}_\infty[h] & \mat{M}_\infty[h]
      \end{array}
      \right]^m.
\end{equation*}
Because $\mat{M}_\infty[h]^m$ satisfies \eqref{eq:infinite-matrix-condition}
for all $m\in\mathbb{N}$,
then also
\begin{equation*}
  \left[
    \begin{array}{cc}
      \mat{P} & \mat{N}\\
      \mat{N} & \mat{P}
    \end{array}
    \right]^m
\end{equation*}
satisfies \eqref{eq:infinite-matrix-condition} for all $m\in\mathbb{N}$.
By Theorem~\ref{thm:m-positive-convergence}, we have convergence for
the diagonal elements of the $2\,m$ powers. By
Lemma~\ref{lemma:convergence-of-norm}, we have
\begin{equation*}
  \mat{e}_i^\top\,
    \left[
    \begin{array}{cc}
      \mat{P}_n & \mat{N}_n\\
      \mat{N}_n & \mat{P}_n
    \end{array}
    \right]^m\rightarrow
  \mat{e}_i^\top\,
    \left[
    \begin{array}{cc}
      \mat{P} & \mat{N}\\
      \mat{N} & \mat{P}
    \end{array}
    \right]^m
\end{equation*}
as $n\rightarrow\infty$.
The convergence for polynomially
bounded functions follows from
Theorem~\ref{thm:dominated-convergence}.
\end{proof}

\subsection{Improper integrals on finite endpoints of an interval}
\label{sec:improper-integrals-on-finite-endpoints-of-an-interval}
In this section, we analyse the convergence for integrals that are
singular on a finite endpoint. A basic requirement for the convergence
is that a singular endpoint of the function is not a node of
the numerical integration rule. It is already known that the nodes of numerical integration, that is, the eigenvalues of
$\mat{M}_n[g]$ lie in the closed interval
$[\operatorname{ess\,inf} g,\operatorname{ess\,sup} g]$
\cite[Theorem~1]{Sarmavuori+Sarkka:2019}.
The \textit{essential infimum}
$\operatorname{ess\,inf} g$ and the
\textit{essential supremum} $\operatorname{ess\,sup} g$ are
the endpoints of the convex hull of the essential range
$\mathcal{R}(g)$.

In the case of an improper integral, where the integrand is singular
at an endpoint of integration, we do not want the endpoint to be
an evaluation point for the numerical integral, that is,
if the outside function $f$ is singular at $\operatorname{ess\,inf} g$
or $\operatorname{ess\,sup} g$, we do not want $\operatorname{ess\,inf} g$
or $\operatorname{ess\,sup} g$ to be an eigenvalue of $\mat{M}_n[g]$.
Fortunately, we have the following theorem.
\begin{theorem}
  \label{Thm:endpoints}
For a measurable real function $g$,
$\operatorname{ess\,inf} g$ (or $\operatorname{ess\,sup} g$)
cannot be an eigenvalue of
the matrix $\boldsymbol{M}_n[g]$ unless
\begin{enumerate}
\item $\operatorname{ess\,inf} g$ (or $\operatorname{ess\,sup} g$) is an eigenvalue of $\mathsf{M}[g]$ and
\item the corresponding eigenfunction is $\phi=\sum_{i=0}^n c_i\,\phi_i$.
\end{enumerate}
\end{theorem}
\begin{proof}
We define a function $\tilde{g}=g-\operatorname{ess\,\inf} g$.
Now $\operatorname{ess\,\inf} \tilde{g}=0$.
If 0 is an eigenvalue of $\boldsymbol{M}_n[\tilde{g}]$ then
for some non-zero $\boldsymbol{c}\in\mathbb{C}^{n+1}$
\begin{equation*}
0=\boldsymbol{c}^*\,\boldsymbol{M}_n[\tilde{g}]\,\boldsymbol{c}.
\end{equation*}
This is possible only if
\begin{equation*}
0=\int_{\Omega}\tilde{g}(\boldsymbol{x})\,
\left|
\sum_{i=0}^n c_i\,\phi_i(\boldsymbol{x})
\right|^2\,w(\boldsymbol{x})\,d\boldsymbol{x},
\end{equation*}
which is possible only if 0 is an eigenvalue of $\mathsf{M}[\tilde{g}]$ and
$\phi=\sum_{i=0}^n c_i\,\phi_i$
is the corresponding eigenfunction. In that case,
$\operatorname{ess\,\inf} g$ is an eigenvalue of $\mathsf{M}[g]$ and
$\phi=\sum_{i=0}^n c_i\,\phi_i$ the corresponding eigenfunction.
The proof for $\operatorname{ess\,sup} g$ is similar with function
$\tilde{g}=\operatorname{ess\,sup} g - g$.
\end{proof}
For example, when $g(x)=x$ and when the basis functions are polynomials, we have the
well known property of Gaussian quadrature that the endpoints are
not nodes of the Gaussian quadrature rule \cite[Theorem~1.46]{Gautschi:2004}.
For a Gaussian quadrature rule,
this property follows from two facts: multiplication operator
$\operator{M}[g]$ does not have
eigenvalues for $g(x)=x$
and a finite order polynomial cannot be a characteristic function
for a set with positive measure.

Once we have shown that the integrand of an improper integral can not
have singularities at the eigenvalues of $\mat{M}_n[g]$, we
can use Jensen's operator inequality for showing the
convergence of the numerical integral. First, we extend the notion of
operator convexity for nonnegative functions of unbounded self-adjoint
infinite matrices.
\begin{definition}
  \label{def:infinite-matrix-convex-function}
  Let $\mathcal{M}$ be the von Neumann algebra of bounded infinite matrices,
  that is, bounded operators on $\ell^2$. A nonnegative function
  $f$ of an unbounded self-adjoint infinite matrix $\mat{A}$ is an unbounded
  operator $f(\mat{A})$ corresponding to an element
  $m_{f(\mat{A})}\in\widehat{\mathcal{M}}_+$.
  A nonnegative function $f$ defined on an interval $I$ is said to be
  \textit{operator convex for infinite matrices}
  if
  for each $\lambda\in[0,1]$, we have
  \begin{equation*}
    f(\lambda\,\mat{A}+(1-\lambda)\,\mat{B}){\scriptstyle\overset{qf}\leq} \lambda\,
    f(\mat{A})+(1-\lambda)\,f(\mat{B}),
  \end{equation*}
  that is, we have
  $m_{f(\lambda\,\mat{A}+(1-\lambda)\,\mat{B})},m_{f(\mat{A})},
  m_{f(\mat{B})}\in\widehat{\mathcal{M}}_+$, corresponding to operators
  $f(\lambda\,\mat{A}+(1-\lambda)\,\mat{B}),f(\mat{A}),
  f(\mat{B})$ respectively such that
  \begin{equation*}
    m_{f(\lambda\,\mat{A}+(1-\lambda)\,\mat{B})}(\varphi)
    \leq \lambda\,m_{f(\mat{A})}(\varphi)+
    (1-\lambda)\,m_{f(\mat{B})}(\varphi),
  \end{equation*}
  for all
  $\varphi\in\mathcal{M}_*^+$.
\end{definition}

Now we can present a form of Jensen's operator inequality of
\cite[Theorem~2.1]{Hansen+Pedersen:2003}
for unbounded
self-adjoint infinite matrices.

\begin{theorem}
  \label{thm:unbounded-jensen}
  Let a nonnegative function $f$ be operator convex for unbounded infinite
  matrices
  on an interval $I$ as in Definition~\ref{def:infinite-matrix-convex-function}
  and $\mat{A}$ be a self-adjoint infinite matrix
  with $\sigma(\mat{A})\subset I$.
  If $\mat{P}$ is an orthogonal projection to a finite dimensional subspace
  of $\ell^2$ such that $\|\mat{A}\,\mat{P}\,\mat{e}_i\|<\infty$ for all
  $\mat{e}_i$,
  then
  \begin{equation*}
  \mat{P}\,f(\mat{P}\,\mat{A}\,\mat{P})\,\mat{P}
  {\scriptstyle\overset{qf}\leq}
  \mat{P}\,f(\mat{A})\,\mat{P}.
  \end{equation*}
\end{theorem}
\begin{proof}
The proof follows similar steps as the proof in \cite[Theorem~2.1]{Hansen+Pedersen:2003}.
Let $\mat{I}_\infty$ be the identity operator in $\ell^2$ and
$\mat{R}=\mat{I}_\infty-\mat{P}$. Because $\mat{R}$ is an orthogonal projection
to an infinite dimensional subspace of $\ell^2$,
we can write $\mat{R}=\mat{V}\,\mat{V}^*$ and
$\mat{V}^*\,\mat{V}=\mat{I}_{\infty}$
\cite[p.~197]{Reed+Simon:1981}.
We have
$\|\mat{A}\,\mat{R}\,\mat{e}_i\|<\infty$
for all $\mat{e}_i$ because
$\|\mat{A}\,(\mat{I}_{\infty}-\mat{P})\,\mat{e}_i\|\leq
\|\mat{A}\,\mat{e}_i\|+\|\mat{A}\,\mat{P}\,\mat{e}_i\|<\infty$.

We define an operator block matrix
$\mat{X} : \vectorspace{D}(\mat{A})\times \ell^2\mapsto
\ell^2\times\ell^2$
\begin{equation*}
  \mat{X}=
  \left[
    \begin{array}{cc}
      \mat{A} & \mat{0}\\
      \mat{0} & s\,\mat{I}_{\infty}
    \end{array}
    \right],
\end{equation*}
where $s$ is a
scalar in $I$.
We also define two unitary
block matrix operators on $\ell^2\times\ell^2$
\begin{equation*}
  \mat{U}_1
  =
  \left[
    \begin{array}{cc}
      -\mat{V} & \mat{P}\\
      \mat{0}  & \mat{V}^*
    \end{array}
    \right],~~
  \mat{U}_2
  =
    \left[
    \begin{array}{cc}
      \mat{V} & \mat{P}\\
      \mat{0} & \mat{V}^*
    \end{array}
    \right].
\end{equation*}

We define two block matrix operators
\begin{align*}
  \mat{Z}
  &=
  \mat{U}_1^*\,\mat{X}\,\mat{U}_1
  =
  \left[
    \begin{array}{cc}
      \mat{V}^*\,\mat{A}\,\mat{V} &
      -\mat{V}^*\,\mat{A}\,\mat{P}\\
      -\mat{P}\,\mat{A}\,\mat{V} &
      \mat{P}\,\mat{A}\,\mat{P}+s\,\mat{R}
    \end{array}
    \right]\\
  \mat{Y}
  &=
  \mat{U}_2^*\,\mat{X}\,\mat{U}_2
    =
  \left[
    \begin{array}{cc}
      \mat{V}^*\,\mat{A}\,\mat{V} &
      \mat{V}^*\,\mat{A}\,\mat{P}\\
      \mat{P}\,\mat{A}\,\mat{V} &
      \mat{P}\,\mat{A}\,\mat{P}+s\,\mat{R}
    \end{array}
    \right].
\end{align*}
The domain of both operators is the same
$\vectorspace{D}(\mat{Z})=\vectorspace{D}(\mat{Y})=
\vectorspace{D}(\mat{A}\,\mat{V})\times
\vectorspace{D}(\mat{A}\,\mat{P})$.
We can find a countable orthonormal basis of vectors where each basis vector is
in the domain of
both operators $\mat{Z}$ and $\mat{Y}$ and both therefore
have an infinite matrix
representation with respect to the same basis.
We construct this orthonormal basis
by constructing a basis of $\vectorspace{D}(\mat{A}\,\mat{V})$
and a basis of $\vectorspace{D}(\mat{A}\,\mat{P})$.
We find that vectors
$\mat{V}^*\,\mat{e}_i\in\vectorspace{D}(\mat{A}\,\mat{V})$ because
$\|\mat{A}\,\mat{V}\,\mat{V}^*\,\mat{e}_i\|=
\|\mat{A}\,\mat{R}\,\mat{e}_i\|<\infty$.
Vectors $\mat{V}^*\,\mat{e}_i$ also span
$\vectorspace{D}(\mat{A}\,\mat{V})$ because each vector
$\mat{u}\in\vectorspace{D}(\mat{A}\,\mat{V})$ is such that
$\mat{V}\,\mat{u}\in\mat{V}\,\ell^2=
\mat{R}\,\ell^2=\mat{V}\,\mat{V}^*\,\ell^2$.
We notice that $\vectorspace{D}(\mat{A}\,\mat{P})$ is spanned by
vectors $\mat{e}_i$
because
$\|\mat{A}\,\mat{P}\,\mat{e}_i\|<\infty$.
Thus, the operators $\mat{Z}$ and $\mat{Y}$ have an infinite matrix
representation
in a countable basis with linearly independent basis vectors
\begin{equation*}
  \phi_{i,j}=
  \left[
    \begin{array}{c}
      \mat{V}^*\,\mat{e}_i\\
      \mat{e}_j
      \end{array}
    \right]
\end{equation*}
that can be orthonormalized to a countable orthonormal basis $\psi_{i,j}$.
Because $\mat{Z}$ and $\mat{Y}$ have an infinite matrix
representation
in the same orthonormal basis, their sum also has an infinite matrix
representation
in the same basis and we have
\begin{equation*}
  \frac{1}{2}\,\mat{Z}+\frac{1}{2}\,\mat{Y}=
  \left[
    \begin{array}{cc}
      \mat{V}^*\,\mat{A}\,\mat{V} & \mat{0}\\
      \mat{0}                     & \mat{P}\,\mat{A}\,\mat{P}+
      s\,\mat{R}
    \end{array}
    \right].
\end{equation*}

Operators $f(\mat{Z})$ and $f(\mat{Y})$ are nonnegative definite
\begin{align*}
  \mat{0}{\scriptstyle\overset{qf}\leq} f(\mat{Z})
  &=
  \mat{U}_1^*\,f(\mat{X})\,\mat{U}_1
  =
  \left[
    \begin{array}{cc}
      \mat{V}^*\,f(\mat{A})\,\mat{V} & -
      \mat{V}\,f(\mat{A})\,\mat{P}\\
      -\mat{P}\,f(\mat{A})\,\mat{V} & \mat{P}\,f(\mat{A})\,\mat{P}+f(s)\,\mat{R}
    \end{array}
    \right]\\
  \mat{0}{\scriptstyle\overset{qf}\leq} f(\mat{Y})
  &=
  \mat{U}_2^*\,f(\mat{X})\,\mat{U}_2
  =
  \left[
    \begin{array}{cc}
      \mat{V}^*\,f(\mat{A})\,\mat{V} & \mat{V}\,f(\mat{A})\,\mat{P}\\
      \mat{P}\,f(\mat{A})\,\mat{V} & \mat{P}\,f(\mat{A})\,\mat{P}+f(s)\,\mat{R}
    \end{array}
    \right].
\end{align*}
We can take the von Neumann algebra $\mathcal{M}$ of all bounded operators
on $\ell^2\times\ell^2$ and there exist
$m_{f(\mat{Z})},m_{f(\mat{Y})},
m_{f\left(\frac{1}{2}\,\mat{Z}+\frac{1}{2}\,\mat{Y}\right)}
\in\widehat{\mathcal{M}}_+$
because operators $\mat{Z},\mat{Y}$ and their sum have infinite matrix
representation
in the same basis. Therefore,
we have
$m_{f\left(\frac{1}{2}\,\mat{Z}+\frac{1}{2}\,\mat{Y}\right)}(\varphi)\leq
\frac{1}{2}\,m_{f(\mat{Z})}(\varphi)+\frac{1}{2}\,m_{f(\mat{Y})}(\varphi)$
for all $\varphi\in\mathcal{M}_*^+$ and also
\begin{align*}
  \left[
    \begin{array}{cc}
      f(\mat{V}^*\,\mat{A}\,\mat{V}) & \mat{0}\\
      \mat{0}                        & f(\mat{P}\,\mat{A}\,\mat{P}+s\,\mat{R})
      \end{array}
      \right]
  &= f
  \left(
  \frac{1}{2}\,\mat{Z}+\frac{1}{2}\,\mat{Y}
  \right)\\
  &
  {\scriptstyle\overset{qf}\leq}
  \frac{1}{2}\,f(\mat{Z})+\frac{1}{2}\,f(\mat{Y})\\
  &=
  \left[
    \begin{array}{cc}
      \mat{V}^*f(\mat{A})\,\mat{V} & \mat{0}\\
      \mat{0} & \mat{P}\,f(\mat{A})\,\mat{P}+f(s)\,\mat{R}
    \end{array}
    \right].
\end{align*}
This inequality can hold only if we have
\begin{align*}
  f(\mat{P}\,\mat{A}\,\mat{P}+s\,\mat{R})
  &{\scriptstyle\overset{qf}\leq}
  \mat{P}\,f(\mat{A})\,\mat{P}+f(s)\,\mat{R}\\
  \mat{P}\,f(\mat{P}\,\mat{A}\,\mat{P})\,\mat{P}+
  f(s)\,\mat{R}
  &{\scriptstyle\overset{qf}\leq}
  \mat{P}\,f(\mat{A})\,\mat{P}+f(s)\,\mat{R}\\
  \mat{P}\,f(\mat{P}\,\mat{A}\,\mat{P})\,\mat{P}
  &{\scriptstyle\overset{qf}\leq}
  \mat{P}\,f(\mat{A})\,\mat{P},
\end{align*}
because $\mat{R}$ and $\mat{P}$ are orthogonal projections to
orthogonal subspaces and we therefore have the identity
$f(\mat{P}\,\mat{A}\,\mat{P} +s\,\mat{R})=
(\mat{P}+\mat{R})\,f(\mat{P}\,\mat{A}\,\mat{P}+s\,\mat{R})\,\mat{P}
+
(\mat{P}+\mat{R})\,f(\mat{P}\,\mat{A}\,\mat{P}+s\,\mat{R})\,\mat{R}
=
\mat{P}\,f(\mat{P}\,\mat{A}\,\mat{P})\,\mat{P}
+
f(s)\,\mat{R}
$.
\end{proof}

We can use the unbounded infinite matrix Jensen's inequality for the following
theorem that establishes monotone convergence of the numerical integral.
\begin{theorem}
  \label{theorem:operator-convex-monotone}
  Let $\operator{A}$ be a self-adjoint operator
  with infinite matrix representation with respect to orthonormal basis vectors
  $\{\phi_i\}_{i=0}^\infty$ and vector $\phi=\sum_{k=0}^m v_k\,\phi_k$.
  Let nonnegative function $f$ be a operator convex for infinite
  matrices according to Definition~\ref{def:infinite-matrix-convex-function}
  on an interval $I$ that
  contains
  $\sigma(\operator{A})$.
  Let finite approximations
  %
  %
  be
  \begin{align*}
    [\mat{A}_n]_{i,j}&=\langle \phi_i,\operator{A}\,\phi_j\rangle,\\
    [\mat{v}_n]_i&=v_i
  \end{align*}
  for $i,j=0,1,\ldots,n$.
  Then we have
 $   \mat{v}_n^*\,f(\mat{A}_n)\,\mat{v}_n
    \rightarrow
    \langle \phi,f(\operator{A})\,\phi\rangle$.
  %
\end{theorem}
\begin{proof}
  We
  define projection operators $\operator{P}_n$ through relation
  \begin{equation*}
    \operator{P}_n\,\phi=\sum_{k=0}^n\langle\phi_k,\phi\rangle\,\phi_k.
  \end{equation*}
  For $n>m$,
  we use Jensen's operator inequality of
  Theorem~\ref{thm:unbounded-jensen} to obtain
    %
  \begin{align*}
    \mat{v}_n^*\,f(\mat{A}_n)\,\mat{v}_n
    &=
    \langle\phi,\operator{P}_n\,f(\operator{P}_n\,\operator{A}\,\operator{P}_n)\,\operator{P}_n\,\phi\rangle
    %
    \leq
    \langle\phi,\operator{P}_n\,
    f(\operator{A})\,\operator{P}_n\,\phi\rangle
    %
    =
    \langle\phi,f(\operator{A})\,\phi\rangle.
  \end{align*}
  We then have the convergence by Lemma~\ref{lemma:bounded-convergence}.
\end{proof}

As an example, we can consider the convergence for functions that are bounded
by an operator convex function $h(x)=a+\frac{b}{x^p}$ where $a,b\geq 0$ and
$0\leq p\leq 1$.
\begin{theorem}
  \label{thm:finite-improper-integral}
  Let a self-adjoint operator $\operator{A}$ have an infinite matrix
  representation
$\mat{A}_\infty$ with respect to orthonormal basis vectors
$\{\phi_i\}_{i=0}^\infty$.
Let a function $f$ have a singularity at endpoint $c$ of the spectrum of
$\mat{A}_\infty$
such that $|f(x)|\leq a + \frac{b}{|x-c|^p}$ where $a,b\geq 0$ and
$0\leq p \leq 1$ and
$\mat{A}_\infty$ has no
eigenvalue and eigenvector pair such that the endpoint $c$ is the eigenvalue
and the corresponding eigenvector is in $\ell^2_0$. Let $f$ also be
Riemann--Stieltjes integrable with respect to the spectral family of
$\operator{A}$ on any finite interval that does not contain $c$.
Let $\mat{u},\mat{v}\in\ell^2_0$ and
vectors $\psi_1=\sum_{i=0}^\infty [\mat{u}]_i\,\phi_i$ and
$\psi_2=\sum_{i=0}^\infty[\mat{v}]_i\,\phi_i$
be such that $\langle \psi_k,(a+b\,|\operator{A}-c|^{-p})\,\psi_k\rangle<\infty$
for both $k=1,2$.
Let finite
approximations be
\begin{align*}
  [\mat{A}_n]_{i,j}&=\langle \phi_i,\operator{A}\,\phi_j\rangle,\\
  [\mat{u}_n]_i &= [\mat{u}]_i,\\
  [\mat{v}_n]_i &= [\mat{v}]_i,
\end{align*}
for $i,j=0,1,\ldots,n$. Then
\begin{equation*}
  \mat{u}_n^*\,f(\mat{A}_n)\,\mat{v}_n\rightarrow
  \langle \psi_1,f(\operator{A})\,\psi_2\rangle.
\end{equation*}
\end{theorem}
\begin{proof}
  It is enough to consider convergence for
  $  \mat{v}_n^*\,f(\mat{A}_n)\,\mat{v}_n\rightarrow
  \langle\psi_2,f(\operator{A})\,\psi_2\rangle$.
  The more general convergence follows
  from polar decomposition
  \begin{align*}
    4\,\langle\phi,f(\operator{A})\,\psi\rangle
    = &
    \langle\phi+\psi,f(\operator{A})\,(\phi+\psi)\rangle
    -\langle\phi-\psi,f(\operator{A})\,(\phi-\psi)\rangle\\
    &
    +i\,\langle\phi+i\,\psi,f(\operator{A})\,(\phi+i\,\psi)\rangle
    -i\,\langle\phi-i\,\psi,f(\operator{A})\,(\phi-i\,\psi)\rangle.
  \end{align*}
  If
  $\phi,\psi\in\operatorname{span}\{\phi_i\}_{i=0}^m=\{\sum_{i=0}^mc_i\,\phi_i\}$
  then also
  $\phi+\psi,\phi-\psi,\phi+i\,\psi,\phi-i\,\psi\in\operatorname{span}\{\phi_i\}_{i=0}^m$.

  It is also enough to show the convergence
  when the lower endpoint of $\sigma(\mat{A}_\infty)$ is 0
    and $\mat{A}_\infty$ has no eigenvalue at 0 such that the corresponding
    eigenvector is in $\ell^2_0$. Then we can select
    $h(x)=\frac{1}{x^p}$ which is operator convex function on $(0,\infty)$.
    The convergence
    for more general $h(x)=a+\frac{b}{|x-c|^p}$ follows by simple
    transformations.
    The convergence for $h(x)$ follows from
    Theorem~\ref{theorem:operator-convex-monotone}.
    The convergence for $|f(x)|\leq h(x)=a+\frac{b}{|x-c|^p}$ follows
    from Theorem~\ref{thm:scalar-dominated-convergence}.
\end{proof}

For an improper numerical integral on a finite endpoint,
we can present the following
theorem.
\begin{theorem}
  Let orthonormal functions $\phi_0=1,\phi_1,\phi_2,\ldots$ be dense
  in $\vectorspace{L}_w^2(\Omega)$. Let a function $g$ satisfy
  \eqref{eq:mult-mat}
  and let matrices $\mat{M}_n[g]$ have
  elements as in \eqref{eq:matrix-elements}.
  Let $c$ be such $\operatorname{ess\,sup} g$ or $\operatorname{ess\,inf} g$
  that it is not an eigenvalue of $\mat{M}_n[g]$, that is, $c$ is not an
  eigenvalue of $\operator{M}[g]$ or if it is, then the corresponding
  eigenfunction is not a finite linear combination of the basis functions
  $\phi_i$.
  Let a function $f$ be Riemann--Stieltjes integrable with respect to
  the spectral family of $\operator{M}[g]$ on any finite interval
  not containing $c$.
  Let $f$ be  bounded so that
  $|f(x)|\leq a+\frac{b}{|x-c|^p}$ for some $a,b\geq 0$ and $0\leq p\leq 1$
  that satisfy
  \begin{equation*}
    \int_\Omega
    \left(
    a+\frac{b}{|g(\mat{x})-c|^p}
    \right)
    \,w(\mat{x})\,d\mat{x}<\infty.
  \end{equation*}
  Then
  \begin{equation*}
    \lim_{n\rightarrow\infty}[f(\mat{M}_n[g])]_{0,0}=
    \int_\Omega f(g(\mat{x}))\,w(\mat{x})\,d\mat{x}.
  \end{equation*}
\end{theorem}
\begin{proof}
This follows from Theorem~\ref{thm:finite-improper-integral}.
\end{proof}

\section{Conclusions}
We have proved several theorems on the convergence of the finite matrix
approximations of multiplication operators
for numerical integration. The theorems also apply generally to
self-adjoint operators.

For integration on unbounded interval, we have shown several theorems
for strong resolvent convergence.
We  have proved strong resolvent convergence of finite matrix truncations
of self-adjoint operators with infinite matrix representation.
We have proved that strong resolvent convergence for sequence of self-adjoint
operators implies strong convergence for Riemann--Stieltjes integrable functions
and strong resolvent convergence for unbounded functions of
the operators.
We have proved that strong resolvent convergence for sequence of
self-adjoint operators implies convergence of quadratically bounded functions
of the operators for certain quadratic forms that can be used in numerical
integration.
We have proved that for multiplication operators of polynomials
or for self-adjoint operators with nonnegative coefficients, the convergence
results can be extended to higher order functions than
quadratically bounded functions.
We have also demonstrated numerically how the convergence can fail
for very fast growing functions if some matrix coefficients are negative.
We have shown a way to improve the convergence by using weights that
are derived from a bounding function of the integrand.

We have also proved theorems for improper integrals with finite
endpoints.
We have shown that endpoints of integration are not eigenvalues
of the finite matrix approximation of the multiplication operator and therefore
not nodes of the matrix method of integration which is also a property
of Gaussian quadrature. We have shown that numerical integrals converge
for improper integrals if the integrand function is bounded by an
operator convex function.

\section*{Acknowledgements}
The work was supported by Academy of Finland.

redpajama_set_name: RedPajamaArXiv