Over-parametrization
The Matlab code you're referencing is using an over-parametrization with constraints instead of an unconstrained parametrization. What @Frank Harrel said about "[methods that use 4 parameters per segment] would allow for discontinuities in the derivatives of the function and a zero-oder discontinuity" is only true for unconstrained parametrizations, so it does not apply to this code. You can still write all twice differentiable piecewise polynomials with four parameters on each interval if you require the parameters to satisfy certain additional constraints. Let's look at an example:
Example of Two Different Representations of a Piecewise Polynomial
Denote the indicator function of a set $T$ by
$$
I_T(x) = \begin{cases}
1, & \text{ if } x \in T \\
0, & \text{ if } x \notin T \\
\end{cases},
$$
and suppose we have two knots $\xi_1$ and $\xi_2$. Let $A = (- \infty, \xi_1]$, $B = (\xi_1, \xi_2]$, and $C = (\xi_2, \infty)$. Then we could define some functions
$$
\begin{align}
g_1 (x) = I_A \ , \ g_2 (x) = x I_A \ & , \ g_3 (x) = x^2 I_A \ , \ g_4 (x) = x^3 I_A \\
g_5 (x) = I_B \ , \ g_6 (x) = x I_B \ & , \ g_7 (x) = x^2 I_B \ , \ g_8 (x) = x^3 I_B \\
g_9 (x) = I_C \ , \ g_{10} (x) = x I_C \ & , \ g_{11} (x) = x^2 I_C \ , \ g_{12} (x) = x^3 I_C \\
\end{align}
$$
Let's also denote the set of all piecewise cubic polynomials with break points at the points $\xi_1$ and $\xi_2$ by $\mathcal{P}$. This set $\mathcal{P}$ is in fact a 12 dimensional vector space with $\{ g_i \}_{i=1}^{12}$ as a basis. If you don't know what all of this means, it essentially means that we need $12$ parameters to describe any piecewise cubic polynomials with break points at the points $\xi_1$ and $\xi_2$, and it means that any such piecewise polynomial $p \in \mathcal{P}$ can be written as $p(x) = \sum_{i=1}^{12} \beta_i g_i (x)$ for 12 parameters $\{ \beta_i \}_{i = 1}^{12}$. If we pick some arbitrary values of the $\beta_i$'s, we can plot one of these functions:
Notice that the function is not continuous (and hence not differentiable, or second differentiable), which makes sense because up to this point we've been talking about the space of all piecewise cubic polynomials with breakpoints $\xi_1$ and $\xi_2$ with no reference to continuity or differentiability. If we want to talk about only the space of continuous piecewise polynomials with breakpoints $\xi_1$ and $\xi_2$, then we can define $\mathcal{P}^0 := \mathcal{P} \cap \mathcal{C}^0 (\mathbb{R})$ where $\mathcal{C}^0 (\mathbb{R})$ means the space of all continuous functions on $\mathbb{R}$.
For a function $p(x) = \sum_{i=1}^{12} \beta_i g_i (x)$ to be continuous, it would need to satisfy two constraints, namely that
$$
p(\xi_1) = \lim_{x \to \xi_1^+} p(x) \quad \text{ and } \quad p(\xi_2) = \lim_{x \to \xi_2^+} p(x)
$$
or in terms of the parameters $\beta_i$:
$$
\begin{align}
\beta_1 + \beta_2 \xi_1 + \beta_3 \xi_1^2 + \beta_4 \xi_1^3 & = \beta_5 + \beta_6 \xi_1 + \beta_7 \xi_1^2 + \beta_8 \xi_1^3 \quad \text{ and }\\
\beta_5 + \beta_6 \xi_2 + \beta_7 \xi_2^2 + \beta_8 \xi_2^3 & = \beta_9 + \beta_{10} \xi_2 + \beta_{11} \xi_2^2 + \beta_{12} \xi_2^3 \\
\end{align}
$$
But $\mathcal{P}^0$ is not a 12-dimensional space like $\mathcal{P}$ is! Essentially, having to satisfy two constraints subtracts two from the dimension to make it a 10-dimensional space. Requiring differentiability would require two more constraints, making $\mathcal{P}^1 := \mathcal{P} \cap \mathcal{C}^1 (\mathbb{R})$ an 8-dimensional space, and requiring twice differentiability would require yet two more constraints, making $\mathcal{P}^2 := \mathcal{P} \cap \mathcal{C}^2 (\mathbb{R})$ a 6-dimensional space. That number should be familiar as the same number of basis functions $\{ h_i \}_{i=1}^6$ you gave in your question, because $\{ h_i \}_{i=1}^6$ is a basis precisely for the 6-dimensional space $\mathcal{P}^2$. Here are the precise constraints in terms of the parameters $\beta_i$ needed to represent $\mathcal{P}^2$:
\begin{align*}
\beta_1 + \beta_2 \xi_1 + \beta_3 \xi_1^2 + \beta_4 \xi_1^3 & = \beta_5 + \beta_6 \xi_1 + \beta_7 \xi_1^2 + \beta_8 \xi_1^3 \\
\beta_5 + \beta_6 \xi_2 + \beta_7 \xi_2^2 + \beta_8 \xi_2^3 & = \beta_9 + \beta_{10} \xi_2 + \beta_{11} \xi_2^2 + \beta_{12} \xi_2^3 \\
\beta_2 + 2 \beta_3 \xi_1 + 3 \beta_4 \xi_1^2 & = \beta_6 + 2 \beta_7 \xi_1 + 3 \beta_8 \xi_1^2 \\
\beta_6 + 2 \beta_7 \xi_2 + 3 \beta_8 \xi_2^2 & = \beta_{10} + 2 \beta_{11} \xi_2 + 3 \beta_{12} \xi_2^2 \\
2 \beta_3 + 6 \beta_4 \xi_1 & = 2 \beta_7 + 6 \beta_8 \xi_1 \\
2 \beta_7 + 6 \beta_8 \xi_2 & = 2 \beta_{11} + 6 \beta_{12} \xi_2 .
\end{align*}
The first two are the continuity constraints from before, the next two are the two differentiability constraints for the points $\xi_1$ and $\xi_2$, and the last two are the second differentiability constraints.
But now we have two different ways of representing functions $p(x) \in \mathcal{P}^2$: we can either write them in terms of the functions $g_i (x)$ as $p(x) = \sum_{i=1}^{12} \beta_i g_i (x)$ where the $\beta_i$ are forced to satisfy the six constraints above, or we can write $p(x) = \sum_{i=1}^{6} \alpha_i h_i (x)$ for some different parameters $\alpha_i$ are aren't constrained at all.
This is the crux of the matter. The Matlab code is doing something similar using the 12 functions $g_i$ and then requiring the parameters to satisfy certain constraints. The set of functions $\{ g_i \}_{i=1}^{12}$ is no longer a basis for this space, because they won't be linearly independent, but you can call this set of functions a spanning set, redundant basis, or over-complete basis (the latter two are not actually bases at all, since they are not linearly independent). This language which is used frequently when talking about frames in signal processing and in dictionary learning.
B-Splines and Numerics
In the actual code, you won't see any functions like $g_i$ or $h_i$. The basis $\{ h_i \}_{i=1}^6$ is called a truncated power basis, and it has some problems. These start getting somewhat technical, so see Carl de Boor, A Practical Guide to Splines Chapter VIII "Example: the truncated power basis can be bad" for details if you are very interested in that. The solution is to use B-splines instead, which are a different basis for the same space with better numerical properties. The code you're referencing uses B-spline in its implementation, so if you want to fully understand it you'll have to learn about them as well.
Terminology
Regression spline describes the concept of taking a fixed set of knots $\{ \xi_i \}_{i=1}^n$ and then fitting them to data according to some criteria, often involving least squares (the definition here is not given on Wikipedia, so see Hastie, Tibshirani, and Friedman The Elements of Statistical Learning Section 5.2 for this definition). Spline for regression and spline fitting also describe the same concept. These aren't popularly used terms, as you've discovered. You might notice that none of these terms have their own Wikipedia article (all three terms redirect to smoothing spline, which has a small note at the end on about regression splines). Smoothing splines are a particular type of regression splines designed to solve the knot-choosing problem, described below.
You said in your question that "This function is distinct from the idea of spline interpolation in that not every data point is used as a knot, only a few." In practice we can choose as many knots as we want and we can set them to be whatever we want. However, too much freedom can be a bad thing when fitting statistical models because of the over-fitting problem. To avoid this, one solution is to limit the number of knots. However, with a small number of knots, the particular choice of which knots to use becomes increasingly important. We now have an additional problem of making this choice. Sometimes theory can help us decide, in the few cases where there is theoretical reason to set a knot in a particular point. However, in many applications there is no good way to choose the knots, and your results will differ depending on the choice of the knots.
Smoothing splines are a related technique that solves this problem, which you can read about on the Wikipedia page.
Sources
Chapter 5 of Hastie, Tibshirani, and Friedman's The Elements of Statistical Learning is an excellent and short source on this. Carl de Boor's A Practical Guide to Splines (if you can find a copy) is the source on splines. Most of the book deals more with splines for interpolation, but chapter XIV is dedicated to smoothing splines.
