Not sure try to fit $\phi(X,f)$ makes sense: how would you define $X$?
For a linear model $X$ is naturally defined as the beta of a linear regression.
Probably you would rather need to go "one layer deeper in the definition of factors", i.e. to use characteristics of the stock (like its market cap, P/E, dividend yield, etc) that you will name $C_1,\ldots,C_N$ and now the goal is to find a neural network (or any machine learning algo) such that
$$R = \sum_{1\leq i\leq F} X_i\cdot\psi_i(C_1,\ldots,C_N)+\epsilon.$$
This is in brief what is proposed in "Deep learning in asset pricing" by Chen, Pelger, and Zhu (2020).
That being said, what you have in mind now is specifically to focus on out of sample risk prediction. This can be done by adding a term in your loss function. What you note now is that in this better formulation, the risk of your portfolio continues to have this kind of shape
$$Var(R_P) = w^T X^T \Omega X w + Var(\epsilon)$$
where the expression of $\Omega$ is now different:
$$\Omega := Cov(\psi_i(C_1,\ldots,C_N),\psi_j(C_1,\ldots,C_N))_{i,j=1,...,F}.$$
Note that the covariance is a scalar product, and hence if your learning algorithms $(\psi_i)_i$ are kernel based, you could use the "kernel trick" (I am not saying that it would be easy... I never seriously thought about doing it).
[EDIT following the comments to address the question of the best way to formulate Non-linear factors]
The linear formulation of factor has the advantage of defining simultaneously the factors and the loadings (here I take a case with more than one factor to express the problem a more generic way, and I explicitly figure the time $t$):
$$R_i(t)=\sum_{i,j}X_i f_j(t) +\epsilon_i(t).$$
It is a well-posed problem in the sense that under quite generic assumptions you can
- use a PCA on returns of a lot of stocks $(R_1,\ldots,R_N)$ to find $K$ factors (see for instance Factors That Fit the Time Series and Cross-Section of Stock Returns by Lettau and Pelger)
- use a linear regression in a second stage to find the loadings.
If you want to replace the factors with a non linear counterpart, my advice is
- work on the residuals of the linear factors, ie on $\epsilon$, because at least it will guarantee you that you are doing something beyond the natural linear factors
- either you use any non-linear PCA like Self Organizing Maps, ie you are in a non supervised mode, and in such a case you have no inputs and the data are again $(R_1,\ldots,R_N)$.
In such a case your loss function should be something like explained variance of the cross section of returns like for a PCA
- either you use a *supervised algorithm (the most famous being neural nets, like perceptrons) and you need exogenous inputs, like characteristics of the stock.
In such a case your loss function will be the variance of the residuals of the explained returns.
It is interesting to notice
- that 1 and 2 respectively correspond to PCA vs characteristics of factors in the linear case (ie there is already 2 approaches that corresponds to 2 different philosophies)
- to really exploit the non linearity of the approach you may decide to go beyond L2 criteria. For instance for case 1 you can try to explain not only the variance of the cross section but its potential skewness or any other non L2 property.