3
$\begingroup$

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?

$\endgroup$
14
  • 3
    $\begingroup$ Avoid using pinv if possible and use the mldivide 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$
    – usεr11852
    Commented Oct 16, 2015 at 17:49
  • $\begingroup$ Ok I added my MATLAB (kinda-pseudo) code. $\endgroup$
    – jeff
    Commented Oct 16, 2015 at 18:03
  • 3
    $\begingroup$ @usεr11852 Is correct. You should almost never explicitly invert a matrix in numerical code, instead favoring to use a linear equation solver. I'm not familiar with matlab, but there should be a function like 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$ Commented Oct 16, 2015 at 18:38
  • $\begingroup$ @MatthewDrury thanks ! How can I do this with MATLAB? $\endgroup$
    – jeff
    Commented Oct 16, 2015 at 18:42
  • 2
    $\begingroup$ It looks like its called linsolve: mathworks.com/help/matlab/ref/linsolve.html $\endgroup$ Commented Oct 16, 2015 at 18:44

1 Answer 1

1
$\begingroup$

Based on what is said in the comments it appears you are solving an under-determined system. That is because the number of samples $N$ is smaller than your number of features $D_1$. Because of this issue one will be always faced with rank deficiency of the design matrix; the thread What is rank deficiency, and how to deal with it? gives more information on the matter.

Notation-wise, as discussed is will be better if one uses standard notation where the number of samples $N$ represent the number of rows rather than the number of columns in the design matrix used. As mentioned the general idea one needs to remember when making this change is that $(AB)^T = B^TA^T$.

As mentioned both by me and @Matthew Drury you should seriously avoid using pinv unless it is explicitly stated in the algorithm you use (so the authors say something like "taking the Moore-Penrose pseudoinverse"). You should use mldivide ( linspace is also an option but it is redundant in your scenario, the Computational Science SE on what differences are between linsolve and mldivide? can offer more clarificaitons on the matter).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.