so the matrix A coefficients are in fact functions and not constants, no ... by that logic any LUT transform can be "claimed" also linear
Functions of R,G,B. They are 'constants' once you decide what powers of R, G, and B you are going to use with a given set of certain R, G, and B data. After that one does not change the matrix A. It is the same methodology used in a 3x3 matrix case, where for the matrix A, the 3 columns are R, G, and B themselves directly. BTW, if you plot those R, G, B, they are typically nonlinear, and not linear at all. Even in a 3x3 setting. In the higher dim-setting I mentioned, one uses nonlinear mappings of these R, G, B. However, that has nothing do to with the transformation being linear of the form y = Ax.
A nonlinear LUT would typically not be able to cast in a linear model of the form y = Ax.