8

我有一个很大的 mxn 矩阵，并且我已经确定了线性相关的列。但是，我想知道 R 中是否有一种方法可以根据线性独立列来编写线性相关列。由于它是一个大矩阵，因此无法根据检查进行。

这是我拥有的矩阵类型的玩具示例。

> mat <- matrix(c(1,1,0,1,0,1,1,0,0,1,1,0,1,1,0,1,0,1,0,1), byrow=TRUE, ncol=5, nrow=4)
> mat
      [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0    1    0
[2,]    1    1    0    0    1
[3,]    1    0    1    1    0
[4,]    1    0    1    0    1

这里很明显 x3 = x1-x2, x5=x1-x4。我想知道是否有一种自动化的方法来获得更大的矩阵。

谢谢！

4

1 回答 1

9

我确信有更好的方法，但我想玩这个。我基本上在开始时检查输入矩阵是否是满列秩，以避免在满秩的情况下进行不必要的计算。之后，我从前两列开始并检查该子矩阵是否具有满列等级，如果是，则检查前列，依此类推。一旦我们找到一些不是全列等级的子矩阵，我就会将该子矩阵中的最后一列回归到前一列，这告诉我们如何构造第一列的线性组合以获得最后一列。

我的功能现在不是很干净，可以做一些额外的检查，但至少这是一个开始。

mat <- matrix(c(1,1,0,1,0,1,1,0,0,1,1,0,1,1,0,1,0,1,0,1), byrow=TRUE, ncol=5, nrow=4)


linfinder <- function(mat){
    # If the matrix is full rank then we're done
    if(qr(mat)$rank == ncol(mat)){
        print("Matrix is of full rank")
        return(invisible(seq(ncol(mat))))
    }
    m <- ncol(mat)
    # cols keeps track of which columns are linearly independent
    cols <- 1
    for(i in seq(2, m)){
        ids <- c(cols, i)
        mymat <- mat[, ids]
        if(qr(mymat)$rank != length(ids)){
            # Regression the column of interest on the previous
            # columns to figure out the relationship
            o <- lm(mat[,i] ~ mat[,cols] + 0)
            # Construct the output message
            start <- paste0("Column_", i, " = ")
            # Which coefs are nonzero
            nz <- !(abs(coef(o)) <= .Machine$double.eps^0.5)
            tmp <- paste("Column", cols[nz], sep = "_")
            vals <- paste(coef(o)[nz], tmp, sep = "*", collapse = " + ")
            message <- paste0(start, vals)
            print(message)
        }else{
            # If the matrix subset was of full rank
            # then the newest column in linearly independent
            # so add it to the cols list
            cols <- ids
        }
    }
    return(invisible(cols))
}

linfinder(mat)

这使

> linfinder(mat)
[1] "Column_3 = 1*Column_1 + -1*Column_2"
[1] "Column_5 = 1*Column_1 + -1*Column_4"
于 2012-10-26T15:19:04.383 回答