雖然 pseudo-inverse 比較 stable 同埋通常比直接計 inverse 快，但並冇比 optimized 過嘅 Gaussian elimination 快，因為 pseudo-inverse 係基於 SVD︰pinv uses the singular value decomposition to form the pseudoinverse of A. Singular values along the diagonal of S that are smaller than tol are treated as zeros, and the representation of A becomes...
可以檢查一下 R 嘅 lm.fit if(method != "qr")
warning(gettextf("method = '%s' is not supported. Using 'qr'", method), domain = NA)
chkDots(...)
z <- .Call(C_Cdqrls, x, y, tol, FALSE)if m <= n:
gufunc = _umath_linalg.lstsq_m
else:
gufunc = _umath_linalg.lstsq_n
所以 linear regression 實際上喺 X full column rank 係用緊 QR factorization。順便講埋 statsmodels 嘅 OLS 係用 MLE (呢個 post 並唔係教書，基礎相關教科書都有提 MLE)X = QR
為免爆 memory 一般都係用 reduced QR 而唔係 full QR. 雖然一般 linear algebra 嘅書會話 QR decomposition 即係 Gram-Schmidt algorithm 嘅結果, 只係 rearrange 咗啲 orthogonal basis 同埋 inner product. 
慎防有人問點解唔用 modified Gram-Schmidt 都係因為 accuracy 問題嚟. 當然 Lapack 入面唔止一種方法 
小妹冇記曬冇睇曬
opts.UT = true;
opts.TRANSA = false;
[Q, R] = qr(X, 0);
b_hat = linsolve(R, Q'*y, opts);b_hat_svd = pinv(X)*y;
當 X 係 10^7 by 5 時用 QR 會快 pinv 約 40%, 當 X 係 10^3 by 500 時用 QR 會快 pinv 約 70%, 但 accuracy 差唔多. 現實你都係鳩禁架啦, 小妹主旨得一個: 記得唔好手多多改咗人地 qr 嘅 default option, 都唔好諗住用咩 gradient descent 去做普通 linear regression 咁弱智, 我間唔中見到有人問點解唔用 gradient descent 做 linear regression. 小妹 delete 咗原本筆記部份 details, 所以會睇落有啲位唔自然. 小妹並冇 cover 曬所有方法, 事實上你睇返 Lapack 入面有樣嘢叫 divide-and-conquer SVD, 同 QR 有得比, 所以 SVD 並冇完全 outdate 嘅
多謝捧場
其實 online linear regression 都係會見到 gradient-based method, 不過呢個 post 只係想講普通 full batch 嘅 linear regression, 所以唔準駁嘴
f <- function (b0, b1, b2, b3, sigma) {
# no d is found in this function env., get d from global env.
return(-sum(log(dnorm(d[["y"]],
mean = b0+b1*d[["x1"]]+b2*d[["x2"]]+b3*d[["x3"]],
sd = sigma))))
}
mle(f,
start = list(b0 = 0.3,
b1 = 0.2,
b2 = -0.4,
b3 = 0.5,
sigma = 1),
method = "BFGS") %>%
coeftest()
又或者唔用mle()
mle_t <- function (data) {
neg_ln_L_t <- function (para, y, x1, x2, x3) {
return(-sum(log(dt(y-para[1]-para[2]*x1-para[3]*x2-para[4]*x3,
df = para[5]))))
}
betas <- optim(par = c(0.3, 0.2, -0.4, 0.5, 1),
fn = neg_ln_L_t,
y = data[["y"]], x1 = data[["x1"]], x2 = data[["x2"]], x3 = data[["x3"]],
method = "BFGS")$par
return(betas)
}
mle_t(data = d)
f <- function (b0, b1, b2, b3, sigma) {
# no d is found in this function env., get d from global env.
return(-sum(log(dnorm(d[["y"]],
mean = b0+b1*d[["x1"]]+b2*d[["x2"]]+b3*d[["x3"]],
sd = sigma))))
}
mle(f,
start = list(b0 = 0.3,
b1 = 0.2,
b2 = -0.4,
b3 = 0.5,
sigma = 1),
method = "BFGS") %>%
coeftest()mle_t <- function (data) {
neg_ln_L_t <- function (para, y, x1, x2, x3) {
return(-sum(log(dt(y-para[1]-para[2]*x1-para[3]*x2-para[4]*x3,
df = para[5]))))
}
betas <- optim(par = c(0.3, 0.2, -0.4, 0.5, 1),
fn = neg_ln_L_t,
y = data[["y"]], x1 = data[["x1"]], x2 = data[["x2"]], x3 = data[["x3"]],
method = "BFGS")$par
return(betas)
}
mle_t(data = d)mle_t <- function (data) {
f <- function (b0, b1, b2, b3, df) {
return(-sum(log(dt(data[["y"]]-b0-b1*data[["x1"]]-b2*data[["x2"]]-b3*data[["x3"]],
df = df))))
}
mle(f,
start = list(b0 = 0.3,
b1 = 0.2,
b2 = -0.4,
b3 = 0.5,
df = 1),
method = "BFGS") %>%
coeftest()
}
mle_t(data = d)
因為小妹覺得 linear regression 值得深入睇吓點運作，可以睇埋數學上同實際上嘅同異
米價就靠大家嚟守護其實OLS estimator 確實係要計(X'X)^{-1}
你只係講緊用咩方法計(X'X)^{-1}
