Inverse Prediction with RANSAC with polynomial features

I am trying to fit a Ransac Regressor to a set of datapoints that looks like this:

import numpy as np
input_size = 100
X = np.random.uniform(size=(input_size, 3))* 10
X = X + np.random.randn(input_size, 3)/1.0
Y = np.zeros_like(X)
Y[:, 0] = 2*X[:, 0] * X[:, 0] + X[:, 1]**2 + 5
Y[:, 1] = 5*X[:,1]**2 + X[:, 2]*X[:, 0]
Y[:, 2] = 2*X[:,0]*X[:,2] + 3*X[:, 2]**2


The shape of both X and Y is 100x3

I want to find the coefficients and intercepts that allow me to predict Y from X. At the same time, I want to be able to predict X from Y (using the same coefficients).

In a simple fitting, I can easily get the coefficients matrix M and the intercept B that allow me to go back between X and Y and vice versa:

from sklearn.linear_model import RANSACRegressor

reg = RANSACRegressor(random_state=42)
reg = reg.fit(X, Y)

B = reg.estimator_.intercept_
M = reg.estimator_.coef_

Y = M*X+B
X = np.linalg.inv(M)*(Y-B)


However, due to the nature of the data, I'm interested in finding a fit that also considered polynomial combinations of X. I first generate a new X matrix consisting of all polynomial combinations of the columns of X, using PolynomialFeatures. X will be fitted with a shape of 100x9, with the form [a, b, c, a^2, ab, ac, b^2, bc, c^2].

from sklearn.preprocessing import PolynomialFeatures
reg = RANSACRegressor(random_state=42)
reg = make_pipeline(PolynomialFeatures(2, include_bias = False), reg)
reg = reg.fit(X, Y)


I can still retrieve the intercept and coefficients of the fitting, and compute the pseudo inverse in order to obtain X from Y.

B = reg.estimator_.intercept_
M = reg.estimator_.coef_

Y = M*X+B
X = np.linalg.pinv(M)*(Y-B)


However, X will now have a shape of 100x9, but I would like to get the degree-1 X values. Is this possible?