| 52 | } |
| 53 | |
| 54 | void Information_root_info_scheme::init_yY () |
| 55 | /* Special Initialisation directly from Information form |
| 56 | * Precondition: |
| 57 | * y,Y |
| 58 | * Postcondition: |
| 59 | * R'*R = Y is PSD |
| 60 | * r = inv(R)'*y |
| 61 | */ |
| 62 | { |
| 63 | // Temporary R Matrix for factorisation |
| 64 | const std::size_t n = x.size(); |
| 65 | LTriMatrix LC(n,n); |
| 66 | // Information Root |
| 67 | Float rcond = LdLfactor (LC, Y); |
| 68 | rclimit.check_PD(rcond, "Initial Y not PD"); |
| 69 | |
| 70 | { // Lower triangular Cholesky factor of LdL' |
| 71 | std::size_t i,j; |
| 72 | for (i = 0; i < n; ++i) |
| 73 | { |
| 74 | using namespace std; // for sqrt |
| 75 | LTriMatrix::value_type sd = LC(i,i); |
| 76 | sd = sqrt(sd); |
| 77 | LC(i,i) = sd; |
| 78 | // Multiply columns by square of non zero diagonal. TODO use column operation |
| 79 | for (j = i+1; j < n; ++j) |
| 80 | { |
| 81 | LC(j,i) *= sd; |
| 82 | } |
| 83 | } |
| 84 | } |
| 85 | R = FM::trans(LC); // R = (L*sqrt(d))' |
| 86 | |
| 87 | UTriMatrix RI(n,n); |
| 88 | RI = R; |
| 89 | bool singular = UTinverse(RI); |
| 90 | assert (!singular); (void)singular; |
| 91 | noalias(r) = prod(FM::trans(RI),y); |
| 92 | } |
| 93 | |
| 94 | void Information_root_scheme::update () |
| 95 | /* Recompute x,X from r,R |