Linearly interpolating in L*a*b* space yields negative RGB values? Is RGB's embedding within L*a*b* non-convex?

Question

I am trying to linearly interpolate between colors in CIELAB space starting out and retrieving colors in a linear RGB space.

My workflow is the following:

1. Take 2 colors in linear RGB, let's say (1,0,0) and (0,0,1).
2. Convert those to XYZ space. I am using this matrix. This should be the matrix for RGB with the E-whitepoint.
3. Convert those to L*a*b* using the formulas on wikipedia. For Xn, Yn, Zn, I am using the corresponding XYZ coordinates of RGB(1,1,1), which is (1,1,1) for this matrix. The coordinates are thus scaled so that the whitepoint has a Y of 1.0.
4. Linearly interpolate component-wise between these two L*a*b* elements.
5. Convert the interpolated points to RGB by reversing the above process (using the inverse transform from wikipedia, as well as the inverse matrix for XYZ->RGB.

The problem is, these intermediate colors are outside my RGB space, even though I obviously start out with 2 colors that are inside it. To be specific: the G channel acquires negative values of about -0.058 when interpolating between RGB(1,0,0) and RGB(0,0,1).

If I use the primaries and whitepoint, i.e. the matrix of linear sRGB instead, the same problem occurs.

Is this normal? It seems plausible to me that a the L*a*b* color space might be such that there's no guarantee that the direct line between two RGB-representable colors is also within that representable region. If it's not normal, where could the error be? If it is normal, how should I deal with it?

My code:

import numpy as np


mat = np.array([[0.49,0.31,0.2],[0.17697,0.81240,0.01063],[0.,0.01,0.99]])

def RGB_to_Lab(R,G,B):
    X,Y,Z = (mat@[[R],[G],[B]])[:,0]
    Xn, Yn, Zn = (mat@[[1],[1],[1]])[:,0]
    f = lambda t: t**(1/3) if t>(6/29)**3 else (1/3)*(29/6)**2*t+4/29
    L = 116*f(Y/Yn) - 16
    a = 500*(f(X/Xn)-f(Y/Yn))
    b = 200*(f(Y/Yn)-f(Z/Zn))
    return L, a, b

def Lab_to_RGB(L,a,b):
    f = lambda t:t**3 if t>6/29 else 3*(6/29)**2*(t-4/29)
    Xn, Yn, Zn = (mat@[[1],[1],[1]])[:,0]
    X = Xn*f((L+16)/116+a/500)
    Y = Yn*f((L+16)/116)
    Z = Zn*f((L+16)/116-b/200)
    R,G,B = (np.linalg.inv(mat)@[[X],[Y], [Z]])[:,0]
    return R,G,B


def steps(n):
    Lab = np.array(RGB_to_Lab(1,0,0))
    Lab2 = np.array(RGB_to_Lab(0,0,1))
    for i in range(n+1):
        r = i/n
        tmp = r*Lab2+(1-r)*Lab
        print(Lab_to_RGB(*tmp))

Answer 1

CIE Lab is bounded by the Visible Spectrum, i.e. Standard Observer sensitivity to light, thus interpolating in that space might generate negative RGB values because practically none of the RGB colourspaces, e.g. sRGB, cover its entirety.

To the question of convexity, RGB spaces are not convex in CIE Lab:

sRGB

ITU-R BT.2020

ACES2065-1