As you know, B-splines smooth, run between the data points.
If you want to interpolate, go through the data points exactly,
using a B-spline function, here's how.
In: data points $f_i$ at the integers $i = 0\ 1\ 2 \dots n$
$\ \ \ \ \ $ a cubic B-spline function $B( x; data_i )$ for $0 \leq x \leq n$
Out: a wrapper function $I( x; data_i )$ that interpolates, $I(i) = f_i$ .
Method:
we know that $B(i) = (f_{i-1} + 4 f_i + f_{i+1}) \,/\, 6, i = 0\ 1\ 2 \dots$
or in matrix form,
$\ b = A f$ with $A$ tridiagonal.
Solve $\ A y = f$
then $\ \ B(y) = B( A^{-1} f ) = f$
i.e. $\ \ \ \ \, I( x; f_i ) \equiv B( x; A^{-1} f )$ interpolates.
Notes:
Solving this tridiagonal system is fast and easy.
(You need boundary conditions such as $f_{-1} = f_0$ or $f_{-1} = 0$ ,
and similarly for $f_{n+1}$.)
A function like $(B + 2 I)\, /\, 3$ is between $B()$ and $I()$:
it smooths less than $B$, but is not as sensitive to noise as $I$ .
This particular blend is called an M-N filter, after
Mitchell and Netravali,
Reconstruction Filters in Computer Graphics,
1988, 8p.
M-N splines also have less overshoot than Interpolating splines; see
what-is-the-maximum-overshoot-of-interpolating-splines-in-d-dimensions.
A fancier way of using B-splines to interpolate, using IIR filters, is given in
M. Unser,
Splines: A perfect fit for signal and image processing ,
1999, 17p .
A Bézier curve with six control points
is defined as
\begin{align}
\mathbf{B_6}(t)
&=
\sum _{i=0}^{6}
{6 \choose i}(1-t)^{6-i}t^{i}\,P_i
\tag{1}\label{1}
,
\end{align}
where $P_i$, $i=0,\dots,6$ are the control points of the spline.
Because of the properties of the convex hull of
the Bezier control points,
to get a visual appearance of the straight line
between the points $A,B$,
one can just set
$P_0=A$, $P_6=B$, and place the other five control
points somewhere on the segment $AB$,
so your choice of
$P_0,P_1,P_2=A$,
$P_4,P_5,P_6=B$,
$P_3=\tfrac12\,(A+B)$ will do for that purpose.
However, to get also the linear expression in \eqref{1},
we need to
expand \eqref{1}, in order to get
a representation as
a polynomial of degree $6$
in the standard form
\begin{align}
a_6t^6+a_5t^5+a_4t^4+a_3t^3+a_2t^2+a_1t+a_0
\tag{2}\label{2}
,
\end{align}
where
\begin{align}
a_0&=P_0
,\\
a_1&= 6\,(P_1-P_0)
,\\
a_2&=15\,(P_0-2\,P_1+P_2)
,\\
a_3&=20\,(-P_0+3\,P_1-3\,P_2+P_3)
,\\
a_4&=15\,(P_0- 4 P_1 + 6 P_2 - 4 P_3+ P_4)
,\\
a_5&= 6\,(-P_0+5\,P_1-10\,P_2+10\,P_3-5\,P_4+P_5)
,\\
a_6&=P_0-6\,P_1+15\,P_2-20\,P_3+15\,P_4-6\,P_5+P_6
.
\end{align}
To get a set of control points $P_i$
such that expression \eqref{2} becomes linear in parameter $t$,
we need to make all coefficients $a_2,\dots,a_6$ zero.
The solution then is just
\begin{align}
P_i&=\tfrac16\,(A\cdot(6-i)+B\cdot i)
,\quad i=0,\dots,6
,
\end{align}
that is, all control points are
evenly distributed along the segment $AB$.
Best Answer
It is curious that you want an exact interpolant for measured data. That's rarely what you really want.
In any case, there are infinitely many ways to interpolate your data. You can set up a polynomial in 8-dimensions with appropriate degree and solve for all the coefficients. Most likely, you will have an under-defined system so you won't be able to solve the coefficients exactly, but rather find an infinite solution set.
However, interpolation is not always useful. The price you pay for going exactly through your data points is that you tend lose control on the behavior of the interpolating surface everywhere between the points. And besides, without a (deterministic) sense of what the function should look like between those points, who's to say that your interpolation makes any sense to begin with?
Furthermore, such interpolants are very sensitive to small perturbations. Suppose you had a very small error in one of your measurements... correcting this error can introduce a whole different interpolating surface!
There are things you should ask yourself about how you expect the curve to behave. Does it go off to infinity as your independent variables get large? Or does it vanish? Or neither?
Do you constrain your first derivatives? Second? Third?
If you cannot answer those questions, you will have a harder time finding the "right" answer. But if all you need is an exact interpolant, then there are methods out there. Polynomial neural networks provide one such approach: set up your eight inputs, one output, some arbitrary number of hidden layers/nodes, and simply continue to solve for the weights until your relative error is within some very small tolerance.
Nevertheless, I caution against it. I can not even think of one situation where that is the "right" thing to do.