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?
pinv
if possible and use themldivide
operator. You do not show the way you try to solve this system. Especially if you use ridge regression, ie. you solve $(X^TX + \lambda I)^{-1} (X^T Y)$ it is even harder to get a singular matrix (you amp all the eigenvalues by $\lambda$ essentially. Please give some reproduction code so we can see your estimation procedure. I strongly suspect that the issue stems from there. $\endgroup$solve(A, y)
which returns the solution(s) to the linear equation $Ax = y$. Always favor this. Here's a nice blog on the subject: johndcook.com/blog/2010/01/19/dont-invert-that-matrix $\endgroup$linsolve
: mathworks.com/help/matlab/ref/linsolve.html $\endgroup$