アルマジロその3. 線形代数編

これまでと同様です. 説明しないので, 気になったらコピペしてみてください.

おまじない

library (Rcpp)
library (inline)
library (RcppArmadillo)
library (Matrix) # for sp_mat 

コレスキー分解

# コレスキー分解
src <- 'arma::mat X = arma::randu(5,5);
        arma::mat Y = X.t() * X;
        arma::mat R = arma::chol(Y);
        arma::mat RtR = R.t()* R;
        Y.print("Y = "); 
        RtR.print("RtR = ");
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

行列式

#行列式
src <- 'arma::mat A = arma::randu(3,3);
        double dtA  = arma::det(A);
        //double ldtA = arma::log_det(1,1,A); //よくわからん
        
        return Rcpp::NumericVector::create(Named("dtA") = dtA);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

固有値

# 固有値問題1
# 複素はなぜか動かない....
# 固有値だけ欲しいときも動かない....
src <- '// 実対称行列
        arma::mat X = arma::randu(5,5);
        arma::mat Y = X.t() * X;
        arma::vec eigval;
        arma::mat eigvec;
        arma::eig_sym (eigval, eigvec, Y);

        return Rcpp::List::create(eigval, 
                                  eigvec,
                                   Y);

'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
rtn <- fun ()


# 固有値問題2
src <- '// 一般の正方行列
        arma::mat X = arma::randu(5,5);
        arma::cx_vec eigval;
        arma::cx_mat eigvec;
        arma::eig_gen (eigval, eigvec, X);

        return Rcpp::List::create(eigval, 
                                  eigvec,
                                   X);

'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
rtn <- fun ()



# 固有値問題3
# 実正方行列固有値分解
src <- '// AU = BUD
        arma::mat A = arma::randu(10,10);
        arma::mat B = arma::randu(10,10);
        arma::cx_vec eigval;
        arma::cx_mat eigvec;
        arma::eig_pair (eigval, eigvec, A, B);

        arma::cx_mat AU  = A * eigvec;
        arma::cx_mat BUD = B * eigvec * diagmat (eigval);

        return Rcpp::List::create(AU,BUD);

'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
rtn <- fun ()

逆行列

# 逆行列1
src <- 'arma::mat A = arma::randu(10, 10);
        arma::mat B = arma::inv(A);
        arma::mat C = A.i();
        return Rcpp::List::create(B, C);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()


# 逆行列2
# 正定値対称行列
src <- 'arma::mat A   = arma::randu(10, 10);
        arma::mat AtA = A.t() * A;
        arma::mat B   = arma::inv_sympd(AtA);
        return Rcpp::List::create(AtA, B);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
rtn <- fun ()

# 逆行列3
# ペンローズの一般
src <- 'arma::mat A   = arma::randu(4, 5);
        arma::mat B   = arma::pinv(A);
        B.print();
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

PLU分解

# PLU分解 (PA = LU)
# Pは行置換
src <- 'arma::mat A   = arma::randu(5, 5);
        arma::mat L, U, P;
        arma::lu(L, U, P, A);
        arma::mat B  = P.t() * L * U;
        return Rcpp::List::create(A, B, P);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

主成分分析

# 主成分分析
src <- 'arma::mat A = arma::randu(5, 4);
        arma::mat coeff, score;
        arma::vec latent, tsquared;
        
        arma::princomp(coeff, score, latent, tsquared, A);
        return Rcpp::List::create(coeff, score, latent, tsquared, A);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

QR分解

# QR分解
# 直交行列Qと行標準型R
# qr_econもある
src <- 'arma::mat X = arma::randu(5,5);
        arma::mat Q, R;
        arma::qr (Q, R, X);
        return Rcpp::List::create(Q * R, X);'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

連立方程式

# 連立方程式
src <- 'arma::mat A = arma::randu(5,5);
        arma::mat B = arma::randu(5,5);
        arma::vec b = arma::randu(5);

        arma::vec x = arma::solve (A, b);
        arma::mat X = arma::solve (A, B);
        
        arma::vec x2;
        arma::solve(x2, A, b);

        return Rcpp::List::create(A, b, B, X, x, x2);

'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

特異値分解

# 特異値分解
# svd_econもある
src <- 'arma::mat X = arma::randu(5,5);
        arma::mat U, V;
        arma::vec s;
        svd (U, s, V, X);


        return Rcpp::List::create(X, U * diagmat(s) * V.t());
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()

シルベスター方程式

# sylvester equation
# solve AX + BX + C = 0, X is unknown
src <- 'arma::mat A = arma::randu(5,5);
        arma::mat B = arma::randu(5,5);
        arma::mat C = arma::randu(5,5);
        
        arma::mat X1 = syl(A, B, C);
        
        arma::mat X2;
        syl(X2, A, B, C);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()