|
|
Ublas : |
Subject: Re: [ublas] Patch/proposal for overloading operator * in ublas
From: K.M.A.Chai_at_[hidden]
Date: 2009-08-19 10:28:33
On Wed, 19 Aug 2009, Jesse Perla wrote:
>
> This looks excellent. Do you have any regression tests/examples/docs on what these patches do
> (besides the obvious ones that replace matlab functions of course)?
>
> Also, were you able to deal with the conversion between multi_array slices to ublas?
>
Hi Jesse,
Thanks for the interest. As I said in the previous post, I'm using it
for learning GP models. I've a package (also in the state of flux) for
this at http://homepages.inf.ed.ac.uk/s9810791/gpml-c++.tbz.
Unfortunately, there's also no documentation associated with this.
I'll like to clarify that the only patches are under the directory
umatlab-patches. Much of the package should be seen as "add-on" instead of
patches.
Unfortunately, it doesn't do the conversion between multi_array slices
to ublas. What it can do is load a cell variable from a mat file into a
multi_arrray of ublas vectors/matrices.
I can, however, provide the following code fragment for doing
regularised linear regression, and returning the MSE on a test set.
-----
typedef ublas::vector<double> dvector;
typedef ublas::symmetric_matrix<double, lower, column_major> dsymmat;
typedef umatlab::chol_representation<double, lower> dcholrep;
try {
1 dsymmat XX(usebindings<lower>(trans(alias(X)) * X));
noalias(XX) += regulariser * identity_matrix<double>(XX.size2());
2 dvector alpha(mldivide( dcholrep(XX), trans(X) * y));
dvector ypredict( usebindings( Xtest * alpha ) );
MSE = mean( pow(ypredict - ytest, 2) );
}
catch( bad_rcond& rcond ) {
cout << "#W Bad rcondition: " << rcond.rcond << endl;
MSE = NAN;
}
-------
Annotation:
1. trans(alias(Xtr)) * Xtr says that we are doing trans(X)*X, where the
two X's are the same. "alias" is a new container introduced into the
package to annotate such instances for the compiler. Cf. the "noalias" of
ublas.
1. usebindings will parse the expression given as its arguments and try to
match a BLAS2/BLAS3 operation to it. I believe in this case "syrk" is
invoked. If UMATLAB_NO_USEBINDINGS is defined, then usebindings will be a
an identity operation, and the blas implementation will be used.
2. Notice that trans(Xtr)*Tautr is a valid operation, i.e., transpose
matrix times vector
2. dcholrep re-represent the matrix as its Cholesky factors. TOMS-865
algorithm (patched) for cholesky decomposition is currently used.
2. mldivide does a matrix left divide operation (see matlab mldivide).
Currently, it simply depatches the operation depending on whether the
matrix representation is QRP, cholesky, symmetrix matrix, triangular
matrix, or simply a square matrix. A vector or a matrix is returned
depending on its second argument.
Kian Ming (Adam)
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.