I'm trying to obtain an estimator $f(x)=y$ where $x \in  \mathbb{R}^{D_1}$ and $y \in  \mathbb{R}^{D_2}$, both are column vectors.
So my training set $X$ and $Y$ are data matrices of size $D_1 \times N$ and $D_2 \times N$, respectively, where $N$ is the number of samples, and $D$'s are the input (feature) and output dimensions.

So I want to learn $\beta$ that gives $\beta x \sim y$ in a least-squares fashion. I was doing this in MATLAB simply by `beta_hat = Y * pinv(X);` and it seems like working without a problem. Though I want to ask, is this correct?

**My question:** 

Now I want to implement this without `pinv` because I want to add regularization to it, so I came up with this solution (this is without regularization) : 
$\hat \beta = Y (X^TX)^{-1}X^T$ is this correct? It also works but MATLAB complains about this : 

    Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  2.565271e-20.
And even crashes sometimes. So I think I'm making a mistake somewhere, but where?

Thanks in advance,

Edit
---
Here is what my MATLAB code looks like : (I know there are non-initialized variables like `N`, but I just cropped them out, they are working as expected) : 

    Ntr = round(N * 0.7); % Assign first 70% of the samples as training set
    Trains = [1:Ntr]; Tests = [Ntr+1:N];
    XData = zeros(FeatureSize, N);
    YData = zeros(OutputSize, N);
    
    for n=1:N
    	% Collect the independent data (into the columns of X)
    	XData(:,n) = getFeature(sample(n));
    	% Collect output variable for Train samples : 
    	if find(Trains==n)
    		YData(:,n) = getLabel(sample(n));
    	end
    end % for each sample
    
    % Learn model: 
    if strcmp(RegressionType, 'ordinary')
    	C = YData(:,Trains) * pinv(XData(:,Trains));
    elseif strcmp(RegressionType, 'ordinary_myImplementation')
    	X = XData(:,Trains);
    	Y = YData(:,Trains);
    	C = Y * inv(X'*X)*X'; % this is where the error happens. Isn't this the same with pinv(X) ?
    elseif strcmp(RegressionType, 'ridge')
    	X = XData(:,Trains);
    	Y = YData(:,Trains);
    	C = Y * inv(X'*X + alpha*eye(Ntr,Ntr)) * X';
    else, error('Unknown regression type');    end

    % Apply model on Test samples : 
    YData(:,Tests) = C * XData(:,Tests);

Edit 2
---
After @Matthew Drury's suggestion, I replaced the line `C = Y * inv(X'*X)*X';` to `C = linsolve(X',Y')';`

    

But now I'm getting this error: 

    Warning: Rank deficient, rank = 17, tol =  2.729816e-12.

Is this normal?