OLS: covariance matrix invertibility problem when rows < columns I have read that in OLS when a number of rows (i.e. observations) is smaller than a number of columns (i.e. variables), the covariance matrix $X^{T}X$ cannot be inverted when parameters are being estimated. In other words the square matrix $X^{T}X$ is singular, indicating linear dependency.  
I have tried to reproduce a non-invertible square matrix $X^{T}X$ but I seem to have a mistake somewhere. I would be appreciative if someone helped me understand where the error is. The reproducible code is added below. 
set.seed(2)
a <- round(matrix(rnorm(n = 6, mean = 0, sd = 3), ncol = 3, nrow = 2), 
digits = 0) # simulate a 3x2 matrix 'a'
a
dim(a)


aTRANS <- t(a) # compute transpose of 'a'
aTRANS
dim(aTRANS)


product <- t(a) %*% a # compute covariance matrix 
product
dim(product)

# Now, the inverse of the covariance matrix is 1/det * adjugate 


# install.packages("RConics") package to compute adjugate/adjoint matrix
library(RConics)

adj <- adjoint(product) # compute adjugate of the covariance matrix 
adj


det <- det(product) # determinant of the covaricance matrix 
det

inverse <- (1/det) * adj 
inverse

If the covariance matrix ($X^{T}X$) cannot be inverted, I was under the impression that R should throw an error about the presence of a singular matrix. However, it seems to 
produce the following output.
            [,1]         [,2]        [,3]
[1,] -1.1259e+15  4.37850e+14 -3.7530e+14
[2,]  4.3785e+14 -1.70275e+14  1.4595e+14
[3,] -3.7530e+14  1.45950e+14 -1.2510e+14

Question: why does R not tell me that the covariance matrix is singular but instead it outputs the above?
 A: I am not familiar with the Rconics package, but the problem seems to operate perfectly well using standard matrix functions in base R.  Using these functions shows that your Gramian matrix is singular, and so it does not invert.  (Note that it still has a pseudo-inverse.)
#Create matrix with more rows than columns
A <- matrix(c(-3, 1, 5, -3, 0, 0), nrow = 2, ncol = 3);

#Create Gramian matrix
G <- t(A) %*% A;

#Find determinant of Gramian matrix
det(G);

[1] 0

#Attempt to invert the Gramian matrix
solve(G);

Error in solve.default(G) : 
  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

#Calculate the Moore-Penrose pseudo-inverse
MASS::ginv(G);
      [,1]  [,2] [,3]
[1,] 2.125 1.125    0
[2,] 1.125 0.625    0
[3,] 0.000 0.000    0

A: I don't know about RConics, but running
solve(product)
produces

Error in solve.default(product) : 
    Lapack routine dgesv: system is exactly singular: U[3,3] = 0

(as it should).
A: RConics did not output a warning or error message because you did not ask RConics to compute directly the inverse. You computed the inverse by yourself using tools provided by RConics. Why should RConics raise an error when you compute the adjoint and the determinant?
I ran your code and got that:
     [,1] [,2] [,3]
[1,]  NaN  NaN  NaN
[2,]  NaN  NaN  NaN
[3,]  NaN  NaN  Inf

Anyway, I suspect that the problem in your case were some round-off errors...
