I have implemented a cubic B-Spline interpolation, not approximation, as follows:

import numpy as np
import math
from geomdl import knotvector

def cox_de_boor( d_, t_, k_, knots_):
    if (d_ == 0):
        if ( knots_[k_] <= t_ <= knots_[k_+1]):
            return 1.0
        return 0.0
    
    denom_l = (knots_[k_+d_] - knots_[k_])
    left = 0.0    
    if (denom_l != 0.0):
        left = ((t_ - knots_[k_]) / denom_l) * cox_de_boor(d_-1, t_, k_, knots_)
    
    denom_r = (knots_[k_+d_+1] - knots_[k_+1])
    right = 0.0
    if (denom_r != 0.0):
        right = ((knots_[k_+d_+1] - t_) / denom_r) * cox_de_boor(d_-1, t_, k_+1, knots_)

    return left + right

def interpolate( d_, P_, n_, ts_, knots_ ):
    A = np.zeros((n_, n_))
    
    for i in range(n_):
        for j in range(n_):
            A[i, j] = cox_de_boor(d_, ts_[i], j, knots_)
    
    control_points = np.linalg.solve(A, P_)
    return control_points

def create_B_spline( d_, P_, t_, knots_):
    sum = Vector()   # just a vector class.
    for i in range( len(P_) ):
        sum += P_[i] * cox_de_boor(d_, t_, i, knots_)
    return sum
    
def B_spline( points_ ):
    d     = 3   # change to 2 for quadratic.
    P     = np.array( points_ )
    n     = len( P )
    ts    = np.linspace( 0.0, 1.0, n )
    knots = knotvector.generate( d, n ) # len = n + d + 1
    
    control_points = interpolate( d, P, n, ts, knots)
    
    crv_pnts = []
    for i in range(10):
        t = float(i) / 9
        crv_pnts.append( create_B_spline(d, control_points, t, knots) )
    return crv_pnts


control_points = [ [float(i), math.sin(i), 0.0] for i in range(4) ]
cps = B_spline( control_points )

Result is OK when interpolating 4 points (control vertices):

Result is NOT OK when interpolating 5 points (control vertices):

Result is OK when interpolating 6 points (control vertices):

and so on...

I noticed two things:

1. The spline does not interpolate properly when the number of control vertices is odd.
2. The spline interpolates properly with any number of vertices when the degree becomes quadratic. So, if you change d = 2, in the B_spline function, the curve will interpolate properly for odd and even number of control vertices.

The cox de boor function is correct and according to the mathematical expression, but with a small alteration on the 2nd conditional expression t[i] <= t **<=** t[i+1] (see my previous SO question here for more details). Also, I used numpy to solve the linear system, which also works as expected. Other than np.linalg.solve, I have tried np.linalg.lstsq but it returns the same results.

I honestly do not know where to attribute this abnormal behaviour. What could cause this issue?

The abnormal behavior described is very interesting and its cause is subtle. Basically, the root cause of this behavior is the Cox-DeBoor function implementation with the <= fix.

In my answer to the OP's previous SO question I give a detailed explanation of why this fix is wrong. In short, this implementation constructs basis functions that are "almost-correct" except for inner-knot values. At the intervals outside the inner knots, the basis functions give correct values - including the t=1 value, which was the reason for the fix in the first place (to prevent a zero-row in the interpolation matrix).

Given this knowledge, we can explain how this mysterious abnormal behavior came to be.

The main things to notice are the arrays ts and knots, which are passed as arguments to the interpolate() function within the function B_spline(). The interpolate() function evaluates the basis functions for all ts. If no inner value of ts is equal to an inner knot then all these values will be correct and the result will be correct. However, if there exists any ts[i]==inner_knot, then the value there will be wrong and will ruin the result.

Now, all that remains is to try to explain the two things noted in the question.

First, notice that the knots are clamped and equally spaced - this is what knotvector.generate( d, n ) does. So for d=3 the knots are: [0,0,0,0, 1/(n-3), 2/(n-3),..., (n-4)/(n-3), 1,1,1,1]. For example, for n=5, the knot vector is [0,0,0,0,0.5,1,1,1,1]. More generally, the inner knots (including 0 and 1) can be computed with the command knots[d:-d] = np.linspace(0, 1, n-d+1).

The ts are computed with the command np.linspace( 0.0, 1.0, n ) (i.e., [0, 1/(n-1), 2/(n-1),..., (n-2)/(n-1), 1]).

"The spline does not interpolate properly when the number of control vertices is odd":

This is because when n is odd and d=3, n-d+1 is also odd and 0.5 must then be one of the knots. By the same argument, 0.5 must also be one of the ts. Therefore, we get a wrong answer!

Interestingly, if n is even and d=3, one can prove that there is no common factor between n-3 and n-1 (two consecutive odd numbers do not have a common factor) and therefore no inner value of ts is equal to an inner knot and the result will be correct! For example, for n=6, ts = [0, 1/5, 2/5, 3/5, 4/5, 1] and knots = [..,0, 1/3, 2/3, 1,..], so no fracture i/3 is equal to a fracture j/5.

"The spline interpolates properly with any number of vertices when the degree becomes quadratic":

This is because when d=2 the inner knots will be np.linspace(0, 1, n-1) (i.e., [0,0,0, 1/(n-2), 2/(n-2),..., (n-3)/(n-2), 1,1,1]) and the ts will be np.linspace(0, 1, n). Since there is no common factor between n-2 and n-1, we won't have i/(n-2) == j/(n-1) so always ts != knots and the result will be correct!

I believe this explains this interesting abnormal interpolation behavior.