Tại sao lm hết bộ nhớ trong khi phép nhân ma trận hoạt động tốt cho các hệ số?

Question

Tôi cố gắng để làm ảnh hưởng bất biến hồi quy tuyến tính với R. Dữ liệu của tôi trông giống như

dte yr id v1 v2 
    . . . . . 
    . . . . . 
    . . . . . 

Sau đó tôi quyết định chỉ đơn giản là thực hiện điều này bằng cách làm cho yr một yếu tố và sử dụng lm:

lm(v1 ~ factor(yr) + v2 - 1, data = df) 

Tuy nhiên, điều này dường như hết bộ nhớ. Tôi có 20 cấp độ trong yếu tố của tôi và df là 14 triệu hàng trong đó có khoảng 2GB để lưu trữ, tôi đang chạy trên một máy tính với 22 GB dành riêng cho quá trình này.

Sau đó tôi quyết định thử mọi thứ theo cách cũ thời: tạo biến giả cho mỗi năm tôi t1 để t20 bằng cách thực hiện:

df$t1 <- 1*(df$yr==1) 
df$t2 <- 1*(df$yr==2) 
df$t3 <- 1*(df$yr==3) 
... 

và chỉ cần tính toán:

solve(crossprod(x), crossprod(x,y)) 

này chạy mà không một vấn đề và tạo ra câu trả lời gần như ngay lập tức.

Tôi đặc biệt tò mò về lm mà làm cho nó hết bộ nhớ khi tôi có thể tính toán các hệ số tốt? Cảm ơn.

Answer 1

Không có câu trả lời nào cho đến nay đã chỉ đúng hướng.

Câu trả lời được chấp nhận bởi @idr đang gây nhầm lẫn giữa lm và summary.lm. lm không tính thống kê chẩn đoán ở tất cả; thay vào đó, summary.lm. Vì vậy, ông đang nói về summary.lm.

@Jake 's câu trả lời là một thực tế về sự ổn định số của yếu tố QR và hệ số LU/Choleksy. Câu trả lời của Aravindakshan mở rộng điều này, bằng cách chỉ ra số lượng các hoạt động điểm nổi phía sau cả hai hoạt động (mặc dù như ông nói, ông không tính vào chi phí cho sản phẩm ma trận chéo máy tính). Nhưng, đừng nhầm lẫn số lượng FLOP với chi phí bộ nhớ. Trên thực tế cả hai phương pháp có cùng một bộ nhớ sử dụng trong LINPACK/LAPACK. Cụ thể, lập luận của ông rằng phương pháp QR tốn nhiều RAM hơn để lưu trữ Q là yếu tố không có thật. Dung lượng được nén như được giải thích trong lm(): What is qraux returned by QR decomposition in LINPACK/LAPACK làm rõ cách tính và lưu trữ hệ số QR. Vấn đề tốc độ của v.v. Chol được nêu chi tiết trong câu trả lời của tôi: Why the built-in lm function is so slow in R? và câu trả lời của tôi trên faster lm cung cấp một thói quen nhỏ lm.chol sử dụng phương pháp Choleksy, nhanh gấp 3 lần phương pháp QR.

@Greg câu trả lời/đề xuất cho biglm là tốt, nhưng nó không trả lời câu hỏi. Vì biglm được đề cập, tôi sẽ chỉ ra rằng QR decomposition differs in lm and biglm. biglm tính toán sự phản ánh của chủ hộ sao cho hệ số R kết quả có đường chéo tích cực. Xem Cholesky factor via QR factorization để biết chi tiết. Lý do mà biglm thực hiện điều này, là kết quả R sẽ giống như yếu tố Cholesky, xem QR decomposition and Choleski decomposition in R để biết thông tin. Ngoài ra, ngoài biglm, bạn có thể sử dụng mgcv. Đọc câu trả lời của tôi: biglm predict unable to allocate a vector of size xx.x MB để biết thêm.

---

Sau một Tóm lại, đó là thời gian để gửi câu trả lời của tôi.

Để phù hợp với một mô hình tuyến tính, lm sẽ

1. tạo ra một khung mô hình;
2. tạo ma trận mô hình;
3. gọi lm.fit để xác thực QR;
4. trả lại kết quả của hệ số QR cũng như khung mô hình trong lmObject.

Bạn đã nói khung dữ liệu đầu vào của mình với 5 cột có giá 2 GB để lưu trữ. Với 20 mức độ yếu tố, ma trận mô hình kết quả có khoảng 25 cột lấy 10 GB dung lượng lưu trữ. Bây giờ chúng ta hãy xem cách sử dụng bộ nhớ tăng lên khi chúng ta gọi lm.

• [môi trường toàn cầu] ban đầu bạn có bộ nhớ 2 GB cho khung dữ liệu;
• [lm envrionment] sau đó nó được sao chép vào khung mô hình, chi phí 2 GB;
• [lm môi trường] thì ma trận mô hình được tạo, chi phí 10 GB;
• [lm.fit môi trường] bản sao ma trận mô hình được tạo sau đó được ghi đè bằng hệ số QR, chi phí 10 GB;
• [lm môi trường] kết quả của lm.fit được trả lại, chi phí 10 GB;
• [môi trường toàn cầu] kết quả của lm.fit được trả lại thêm bởi lm, chi phí 10 GB khác;
• [môi trường toàn cầu] khung mô hình được trả về lm, chi phí 2 GB.

Vì vậy, tổng số RAM 46 GB là cần thiết, lớn hơn nhiều so với RAM 22 GB có sẵn của bạn.

Thực tế nếu lm.fit có thể được "gạch chân" thành lm, chúng tôi có thể tiết kiệm 20 GB chi phí. Nhưng không có cách nào để inline một hàm R trong hàm R khác.

Có lẽ chúng ta có thể lấy một ví dụ nhỏ để xem những gì xảy ra xung quanh lm.fit:

X <- matrix(rnorm(30), 10, 3) # a `10 * 3` model matrix 
y <- rnorm(10) ## response vector 

tracemem(X) 
# [1] "<0xa5e5ed0>" 

qrfit <- lm.fit(X, y) 
# tracemem[0xa5e5ed0 -> 0xa1fba88]: lm.fit 

Vì vậy, trên thực tế, X được sao chép khi truyền vào lm.fit. Chúng ta hãy có một cái nhìn vào những gì qrfit có

str(qrfit) 
#List of 8 
# $ coefficients : Named num [1:3] 0.164 0.716 -0.912 
# ..- attr(*, "names")= chr [1:3] "x1" "x2" "x3" 
# $ residuals : num [1:10] 0.4 -0.251 0.8 -0.966 -0.186 ... 
# $ effects  : Named num [1:10] -1.172 0.169 1.421 -1.307 -0.432 ... 
# ..- attr(*, "names")= chr [1:10] "x1" "x2" "x3" "" ... 
# $ rank   : int 3 
# $ fitted.values: num [1:10] -0.466 -0.449 -0.262 -1.236 0.578 ... 
# $ assign  : NULL 
# $ qr   :List of 5 
# ..$ qr : num [1:10, 1:3] -1.838 -0.23 0.204 -0.199 0.647 ... 
# ..$ qraux: num [1:3] 1.13 1.12 1.4 
# ..$ pivot: int [1:3] 1 2 3 
# ..$ tol : num 1e-07 
# ..$ rank : int 3 
# ..- attr(*, "class")= chr "qr" 
# $ df.residual : int 7 

Lưu ý rằng QR ma trận nhỏ gọn qrfit$qr$qr là lớn như mô hình ma trận X. Nó được tạo ra bên trong lm.fit, nhưng trên lối ra của lm.fit, nó được sao chép. Vì vậy, tổng cộng, chúng tôi sẽ có 3 "bản sao" của X:

• nguyên bản trong môi trường toàn cầu;
• ảnh được sao chép vào lm.fit, ghi đè bằng hệ số QR;
• số được trả về lm.fit.

Trong trường hợp của bạn, X là 10 GB, vì vậy chi phí bộ nhớ được kết hợp với riêng lm.fit đã là 30 GB. Hãy để một mình chi phí khác liên quan đến lm.

---

Mặt khác, chúng ta hãy có một cái nhìn tại

solve(crossprod(X), crossprod(X,y)) 

X mất 10 GB, nhưng crossprod(X) chỉ là một ma trận 25 * 25, và crossprod(X,y) chỉ là một vector chiều dài-25. Chúng quá nhỏ so với X, do đó việc sử dụng bộ nhớ không tăng chút nào.

Có thể bạn đang lo lắng rằng một bản sao cục bộ của X sẽ được thực hiện khi gọi crossprod? Không có gì! Không giống như lm.fit thực hiện cả việc đọc và ghi thành X, crossprod chỉ đọc X, do đó không sao chép được thực hiện. Chúng tôi có thể xác minh điều này với ma trận đồ chơi của chúng tôi X bởi:

tracemem(X) 
crossprod(X) 

Bạn sẽ không thấy thông báo sao chép!

---

Nếu bạn muốn có một bản tóm tắt ngắn cho tất cả ở trên, ở đây là:

• chi phí bộ nhớ cho lm.fit(X, y) (hoặc thậm chí .lm.fit(X, y)) là ba lần lớn như rằng cho solve(crossprod(X), crossprod(X,y));
• Tùy thuộc vào ma trận mô hình lớn hơn khung mô hình, chi phí bộ nhớ cho lm lớn gấp 3 ~ 6 lần đối với solve(crossprod(X), crossprod(X,y)). Giới hạn dưới 3 không bao giờ đạt được, trong khi ngưỡng trên 6 đạt được khi ma trận mô hình giống như khung mô hình. Đây là trường hợp khi không có biến yếu tố hoặc các điều khoản "yếu tố-như nhau", như bs() và poly() vv

Answer 2

lm hoạt động nhiều hơn là chỉ tìm các hệ số cho các tính năng nhập của bạn. Ví dụ, nó cung cấp số liệu thống kê chẩn đoán cho bạn biết thêm về các hệ số trên các biến độc lập của bạn bao gồm cả lỗi tiêu chuẩn và giá trị t của mỗi biến độc lập của bạn.

Tôi nghĩ rằng việc hiểu các thống kê chẩn đoán này rất quan trọng khi chạy hồi quy để hiểu mức độ hồi quy hợp lệ của bạn.

Các tính toán bổ sung này gây ra lm chậm hơn đơn giản là giải quyết các phương trình ma trận cho hồi quy.

Ví dụ, bằng cách sử dụng mtcars bộ dữ liệu:

>data(mtcars) 
>lm_cars <- lm(mpg~., data=mtcars) 
>summary(lm_cars) 

Call:               
lm(formula = mpg ~ ., data = mtcars)       

Residuals:              
    Min  1Q Median  3Q  Max      
-3.4506 -1.6044 -0.1196 1.2193 4.6271      

Coefficients:             
      Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.30337 18.71788 0.657 0.5181    
cyl   -0.11144 1.04502 -0.107 0.9161    
disp   0.01334 0.01786 0.747 0.4635    
hp   -0.02148 0.02177 -0.987 0.3350    
drat   0.78711 1.63537 0.481 0.6353    
wt   -3.71530 1.89441 -1.961 0.0633 .    
qsec   0.82104 0.73084 1.123 0.2739    
vs   0.31776 2.10451 0.151 0.8814    
am   2.52023 2.05665 1.225 0.2340    
gear   0.65541 1.49326 0.439 0.6652    
carb  -0.19942 0.82875 -0.241 0.8122    
---               
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2.65 on 21 degrees of freedom   
Multiple R-squared: 0.869,  Adjusted R-squared: 0.8066  
F-statistic: 13.93 on 10 and 21 DF, p-value: 3.793e-07  

Answer 3

Ngoài những gì Idris cho biết, nó cũng có giá trị chỉ ra rằng lm() không giải quyết cho các thông số bằng cách sử dụng phương trình bình thường như bạn minh họa trong của bạn câu hỏi, nhưng thay vì sử dụng phân tách QR, hiệu quả kém hơn nhưng có xu hướng tạo ra nhiều giải pháp chính xác hơn về số lượng.

Answer 4

Bạn có thể muốn xem xét sử dụng gói biglm. Nó phù hợp với các mô hình lm bằng cách sử dụng các khối dữ liệu nhỏ hơn.

Answer 5

Để xây dựng trên quan điểm của Jake. Giả sử hồi qui của bạn đang cố gắng giải quyết: y = Ax (A là ma trận thiết kế). Với quan sát m và n biến độc lập A là ma trận mxn. Sau đó, chi phí của QR là ~ m*n^2. Trong trường hợp của bạn, nó trông giống như m = 14x10^6 và n = 20. Vì vậy, m*n^2 = 14*10^6*400 là một chi phí đáng kể.

Tuy nhiên với phương trình bình thường bạn đang cố đảo ngược A'A ('chỉ thị chuyển vị), hình vuông và có kích thước nxn. Giải quyết thường được thực hiện bằng LU có chi phí n^3 = 8000. Điều này nhỏ hơn nhiều so với chi phí tính toán của QR. Tất nhiên điều này không bao gồm chi phí của ma trận nhân lên.

Hơn nữa nếu quy trình QR cố gắng lưu trữ ma trận Q có kích thước mxm=14^2*10^12 (!), Thì bộ nhớ của bạn sẽ không đủ. Có thể viết QR để không gặp vấn đề này. Thật thú vị khi biết phiên bản QR nào thực sự đang được sử dụng. Và tại sao chính xác cuộc gọi lm hết bộ nhớ.