Testing for linear dependence among the columns of a matrix I have a correlation matrix of security returns whose determinant is zero. (This is a bit surprising since the sample correlation matrix and the corresponding covariance matrix should theoretically be positive definite.)
My hypothesis is that at least one security is linearly dependent on other securities. Is there a function in R that sequentially tests each column a matrix for linear dependence?
For example, one approach would be to build up a correlation matrix one security at a time and calculate the determinant at each step. When the determinant = 0 then stop as you have identified the security who is a linear combination of other securities.
Any other techniques to identify linear dependence in such a matrix are appreciated.
 A: What I'd try to do here for diagnostic purposes is to take the $502\times 480$ matrix (that is, the transpose) and determine the singular values of the matrix (for diagnostic purposes, you don't need the full singular value decomposition... yet). Once you have the 480 singular values, check how many of those are "small" (a usual criterion is that a singular value is "small" if it is less than the largest singular value times the machine precision). If there are any "small" singular values, then yes, you have linear dependence.
A: Here's a straightforward approach: compute the rank of the matrix that results from removing each of the columns.  The columns which, when removed, result in the highest rank are the linearly dependent ones (since removing those does not decrease rank, while removing a linearly independent column does).
In R:
rankifremoved <- sapply(1:ncol(your.matrix), function (x) qr(your.matrix[,-x])$rank)
which(rankifremoved == max(rankifremoved))

A: Not that the answer @Whuber gave really needs to be expanded on but I thought I'd provide a brief description of the math.
If the linear combination $\mathbf{X'Xv}=\mathbf{0}$ for $\mathbf{v}\neq\mathbf{0}$ then $\mathbf{v}$ is an eigenvector of $\mathbf{X'X}$ associated with eigenvalue $\lambda=0$. The eigenvectors and eigenvalues of $\mathbf{X'X}$ are also the eigenvectors and eigenvalues of $\mathbf{X}$, so the eigenvectors of $\mathbf{X'X}$ associated with eigenvalues near $\lambda=0$ represent the coefficients for approximate linear relationships among the regressors. Principal Component Analysis outputs the eigenvectors and eigenvalues of $\mathbf{X'X}$, so you can use the eigenvectors $\mathbf{v}$ associated with small $\lambda$ to determine if linear relationships exist among some of your regressors.
One method of determining if an eigenvalue is appropriately small to constitute collinearity is to use the Condition Indices:
$$
\mathbf{\kappa_j}=\frac{\lambda_{max}}{\lambda_j}
$$
which measures the size of the smallest eigenvalues relative to the largest. A general rule of thumb is that modest multicollinearity is associated with a condition index between 100 and 1,000 while severe multicollinearity is associated with a condition index above 1,000 (Montgomery, 2009). 
It's important to use an appropriate method for determining if an eigenvalue is small because it's not the absolute size of the eigenvalues, it's the relative size of the condition index that's important, as can be seen in an example. Consider the matrix
$$
\mathbf{X'X}=\left[\begin{array}{rr}
0.001 & 0 & 0 \\
0 & 0.001 & 0 \\
0 & 0 & 0.001 \\
\end{array}
\right].
$$
The eigenvalues for this matrix are $\lambda_1=\lambda_2=\lambda_3=0.001$. Although these eigenvalues appear small the condition index is
$$
\mathbf{\kappa}=\frac{\lambda_{max}}{\lambda_{min}}=1
$$
indicating absence of multicolinearity and , in fact, the columns of this matrix are linearly independent.
Citations
Montgomery, D. (2012). Introduction to Linear Regression Analysis, 5th Edition. John Wiley & Sons Inc.
A: I ran into this issue roughly two weeks ago and decided that I needed to revisit it because when dealing with massive data sets, it is impossible to do these things manually. 
I created a for() loop that calculates the rank of the matrix one column at a time. So for the first iteration, the rank will be 1. The second, 2. This occurs until the rank becomes LESS than the column number you are using. 
Very straightforward:
for (i in 1:47) {

  print(qr(data.frame[1:i])$rank) 
  print(i) 
  print(colnames(data.frame)[i])
  print("###") 
}

for() loop breakdown


*

*calculates the rank for the ith column 

*prints the iteration number

*prints the column name for reference 

*divides the console with "###" so that you can easily scroll through


I am sure that you can add an if statement, I don't need it yet because I am only dealing with 50ish columns. 
Hope this helps!
A: The question asks about "identifying underlying [linear] relationships" among variables.
The quick and easy way to detect relationships is to regress any other variable (use a constant, even) against those variables using your favorite software: any good regression procedure will detect and diagnose collinearity.  (You will not even bother to look at the regression results: we're just relying on a useful side-effect of setting up and analyzing the regression matrix.)
Assuming collinearity is detected, though, what next?  Principal Components Analysis (PCA) is exactly what is needed: its smallest components correspond to near-linear relations.  These relations can be read directly off the "loadings," which are linear combinations of the original variables.  Small loadings (that is, those associated with small eigenvalues) correspond to near-collinearities.  An eigenvalue of $0$ would correspond to a perfect linear relation.  Slightly larger eigenvalues that are still much smaller than the largest would correspond to approximate linear relations.
(There is an art and quite a lot of literature associated with identifying what a "small" loading is.  For modeling a dependent variable, I would suggest including it within the independent variables in the PCA in order to identify the components--regardless of their sizes--in which the dependent variable plays an important role.  From this point of view, "small" means much smaller than any such component.)

Let's look at some examples.  (These use R for the calculations and plotting.)  Begin with a function to perform PCA, look for small components, plot them, and return the linear relations among them. 
pca <- function(x, threshold, ...) {
  fit <- princomp(x)
  #
  # Compute the relations among "small" components.
  #
  if(missing(threshold)) threshold <- max(fit$sdev) / ncol(x)
  i <- which(fit$sdev < threshold)
  relations <- fit$loadings[, i, drop=FALSE]
  relations <- round(t(t(relations) / apply(relations, 2, max)), digits=2)
  #
  # Plot the loadings, highlighting those for the small components.
  #
  matplot(x, pch=1, cex=.8, col="Gray", xlab="Observation", ylab="Value", ...)
  suppressWarnings(matplot(x %*% relations, pch=19, col="#e0404080", add=TRUE))

  return(t(relations))
}

Let's apply this to some random data.  These are built on four variables (the $B,C,D,$ and $E$ of the question).  Here is a little function to compute $A$ as a given linear combination of the others.  It then adds i.i.d. Normally-distributed values to all five variables (to see how well the procedure performs when multicollinearity is only approximate and not exact).
process <- function(z, beta, sd, ...) {
  x <- z %*% beta; colnames(x) <- "A"
  pca(cbind(x, z + rnorm(length(x), sd=sd)), ...)
}

We're all set to go: it remains only to generate $B, \ldots, E$ and apply these procedures.  I use the two scenarios described in the question: $A=B+C+D+E$ (plus some error in each) and $A=B+(C+D)/2+E$ (plus some error in each).  First, however, note that PCA is almost always applied to centered data, so these simulated data are centered (but not otherwise rescaled) using sweep.
n.obs <- 80 # Number of cases
n.vars <- 4 # Number of independent variables
set.seed(17)
z <- matrix(rnorm(n.obs*(n.vars)), ncol=n.vars)
z.mean <- apply(z, 2, mean)
z <- sweep(z, 2, z.mean)
colnames(z) <- c("B","C","D","E") # Optional; modify to match `n.vars` in length

Here we go with two scenarios and three levels of error applied to each.  The original variables $B, \ldots, E$ are retained throughout without change: only $A$ and the error terms vary.

The output associated with the upper left panel was
       A  B  C  D  E
Comp.5 1 -1 -1 -1 -1

This says that the row of red dots--which is constantly at $0$, demonstrating a perfect multicollinearity--consists of the combination $0 \approx A -B-C-D-E$: exactly what was specified.
The output for the upper middle panel was
       A     B     C     D     E
Comp.5 1 -0.95 -1.03 -0.98 -1.02

The coefficients are still close to what we expected, but they are not quite the same due to the error introduced.  It thickened the four-dimensional hyperplane within the five-dimensional space implied by $(A,B,C,D,E)$ and that tilted the estimated direction just a little.  With more error, the thickening becomes comparable to the original spread of the points, making the hyperplane almost impossible to estimate.  Now (in the upper right panel) the coefficients are
       A     B     C     D     E
Comp.5 1 -1.33 -0.77 -0.74 -1.07

They have changed quite a bit but still reflect the basic underlying relationship $A' = B' + C' + D' + E'$ where the primes denote the values with the (unknown) error removed.
The bottom row is interpreted the same way and its output similarly reflects the coefficients $1, 1/2, 1/2, 1$.
In practice, it is often not the case that one variable is singled out as an obvious combination of the others: all coefficients may be of comparable sizes and of varying signs.  Moreover, when there is more than one dimension of relations, there is no unique way to specify them: further analysis (such as row reduction) is needed to identify a useful basis for those relations. That's how the world works: all you can say is that these particular combinations that are output by PCA correspond to almost no variation in the data.  To cope with this, some people use the largest ("principal") components directly as the independent variables in the regression or the subsequent analysis, whatever form it might take.  If you do this, do not forget first to remove the dependent variable from the set of variables and redo the PCA!

Here is the code to reproduce this figure:
par(mfrow=c(2,3))
beta <- c(1,1,1,1) # Also can be a matrix with `n.obs` rows: try it!
process(z, beta, sd=0, main="A=B+C+D+E; No error")
process(z, beta, sd=1/10, main="A=B+C+D+E; Small error")
process(z, beta, sd=1/3, threshold=2/3, main="A=B+C+D+E; Large error")

beta <- c(1,1/2,1/2,1)
process(z, beta, sd=0, main="A=B+(C+D)/2+E; No error")
process(z, beta, sd=1/10, main="A=B+(C+D)/2+E; Small error")
process(z, beta, sd=1/3, threshold=2/3, main="A=B+(C+D)/2+E; Large error")

(I had to fiddle with the threshold in the large-error cases in order to display just a single component: that's the reason for supplying this value as a parameter to process.)

User ttnphns has kindly directed our attention to a closely related thread.  One of its answers (by J.M.) suggests the approach described here.
A: Rank, r of a matrix = number of linearly independent columns (or rows) of a matrix. For a n by n matrix A, rank(A) = n => all columns (or rows) are linearly independent.
