  Let $\xi = \{x^j\}_{j=1}^n$ be a grid of $n$ points in the $d$-cube
${\II}^d:=[0,1]^d$, and $\Phi = \{\phi_j\}_{j =1}^n$ a family of $n$ functions
on ${\II}^d$. We define the linear sampling algorithm $L_n(\Phi,\xi,\cdot)$ for
an approximate recovery of a continuous function $f$ on ${\II}^d$ from the
sampled values $f(x^1), ..., f(x^n)$, by $$L_n(\Phi,\xi,f)\ := \ \sum_{j=1}^n
f(x^j)\phi_j$$.
  For the Besov class $B^\alpha_{p,\theta}$ of mixed smoothness $\alpha$
(defined as the unit ball of the Besov space $\MB$), to study optimality of
$L_n(\Phi,\xi,\cdot)$ in $L_q({\II}^d)$ we use the quantity
$$r_n(B^\alpha_{p,\theta})_q \ := \ \inf_{H,\xi} \ \sup_{f \in
B^\alpha_{p,\theta}} \, \|f - L_n(\Phi,xi,f)\|_q$$, where the infimum is taken
over all grids $\xi = \{x^j\}_{j=1}^n$ and all families $\Phi =
\{\phi_j\}_{j=1}^n$ in $L_q({\II}^d)$. We explicitly constructed linear
sampling algorithms $L_n(\Phi,\xi,\cdot)$ on the grid $\xi = \ G^d(m):=
\{(2^{-k_1}s_1,...,2^{-k_d}s_d) \in \II^d : \ k_1 + ... + k_d \le m\}$, with
$\Phi$ a family of linear combinations of mixed B-splines which are mixed
tensor products of either integer or half integer translated dilations of the
centered B-spline of order $r$. The grid $G^d(m)$ is of the size $2^m m^{d-1}$
and sparse in comparing with the generating dyadic coordinate cube grid of the
size $2^{dm}$. For various $0<p,q,\theta \le \infty$ and $1/p < \alpha < r$, we
proved upper bounds for the worst case error $ \sup_{f \in B^\alpha_{p,\theta}}
\, \|f - L_n(\Phi,\xi,f)\|_q$ which coincide with the asymptotic order of
$r_n(B^\alpha_{p,\theta})_q$ in some cases. A key role in constructing these
linear sampling algorithms, plays a quasi-interpolant representation of
functions $f \in B^\alpha_{p,\theta}$ by mixed B-spline series.
