Patterns in static

Apophenia

Functions
Singular value decompositions, determinants, et cetera.

Functions

double apop_det_and_inv (const gsl_matrix *in, gsl_matrix **out, int calc_det, int calc_inv)
 
gsl_matrix * apop_matrix_inverse (const gsl_matrix *in)
 
double apop_matrix_determinant (const gsl_matrix *in)
 
apop_dataapop_matrix_pca (gsl_matrix *data, int const dimensions_we_want)
 
apop_dataapop_dot (const apop_data *d1, const apop_data *d2, char form1, char form2)
 
gsl_vector * apop_numerical_gradient (apop_data *data, apop_model *model, double delta)
 

Detailed Description

This page describes some standard bits of linear algebra that Apophenia facilitates.

See also the printing functions, Assorted printing functions, and the Convenience functions.

Function Documentation

double apop_det_and_inv ( const gsl_matrix *  in,
gsl_matrix **  out,
int  calc_det,
int  calc_inv 
)

Calculate the determinant of a matrix, its inverse, or both, via LU decomposition. The in matrix is not destroyed in the process.

See also
apop_matrix_determinant, apop_matrix_inverse
Parameters
inThe matrix to be inverted/determined.
outIf you want an inverse, this is where to place the matrix to be filled with the inverse. Will be allocated by the function.
calc_det0: Do not calculate the determinant.
1: Do.
calc_inv0: Do not calculate the inverse.
1: Do.
Returns
If calc_det == 1, then return the determinant. Otherwise, just returns zero. If calc_inv!=0, then *out is pointed to the matrix inverse. In case of difficulty, I will set *out=NULL and return NaN.
apop_data* apop_dot ( const apop_data d1,
const apop_data d2,
char  form1,
char  form2 
)

A convenience function for dot products, which requires less prep and typing than the gsl_cblas_dgexx functions.

Second, it makes some use of the semi-overloading of the apop_data structure. d1 may be a vector or a matrix, and the same for d2, so this function can do vector dot matrix, matrix dot matrix, and so on. If d1 includes both a vector and a matrix, then later parameters will indicate which to use.

Parameters
d1the left part of $ d1 \cdot d2$
d2the right part of $ d1 \cdot d2$
form1't' or 'p': transpose or prime d1->matrix, or, if d1->matrix is NULL, read d1->vector as a row vector.
'n' or 0: no transpose. (the default)
'v': ignore the matrix and use the vector.
form2As above, with d2.
Returns
an apop_data set. If two matrices come in, the vector element is NULL and the matrix has the dot product; if either or both are vectors, the vector has the output and the matrix is NULL.
Exceptions
out->error='a'Allocation error.
out->error='d'dimension-matching error.
out->error='m'GSL math error.
NULLIf you ask me to take the dot product of NULL, I return NULL. [May some day change.]
  • Some systems auto-transpose non-conforming matrices. You input a $3 \times 5$ and a $3 \times 5$ matrix, and the system assumes that you meant to transpose the second, producing a $3 \times 5 \cdot 5 \times 3 \rightarrow 3 \times 3$ output. Apophenia does not do this. First, it's ambiguous whether the output should be $3 \times 3$ or $5 \times 5$. Second, your next run might have three observations, and two $3 \times 3$ matrices don't require transposition; auto-transposition thus creates situations where bugs can pop up on only some iterations of a loop.
  • For a vector $cdot$ a matrix, the vector is always treated as a row vector, meaning that a $3\times 1$ dot a $3\times 4$ matrix is correct, and produces a $1 \times 4$ vector. For a vector $cdot$ a matrix, the vector is always treated as a column vector. Requests for transposition are ignored.
  • As a corrollary to the above rule, a vector dot a vector always produces a scalar, which will be put in the zeroth element of the output vector; see the example.
  • If you want to multiply an $N \times 1$ vector $\cdot$ a $1 \times N$ matrix produce an $N \times N$ matrix, then use apop_vector_to_matrix to turn your vectors into matrices; see the example.
  • A note for readers of Modeling with Data: the awkward instructions on using this function on p 130 are now obsolete, thanks to the designated initializer syntax for function calls. Notably, in the case where d1 is a vector and d2 a matrix, then apop_dot(d1,d2,'t') won't work, because 't' now refers to d1. Instead use apop_dot(d1,d2,.form2='t') or apop_dot(d1,d2,0, 't')

Sample code:

/* A demonstration of dot products and various useful
transformations among types. */
#include <apop.h>
double eps=1e-3;//slow to converge series-->large tolerance.
#define Diff(L, R) Apop_assert(fabs((L)-(R)<(eps)), "%g is too different from %g (abitrary limit=%g).", (double)(L), (double)(R), eps);
int main(){
int len = 3000;
gsl_vector *v = gsl_vector_alloc(len);
for (double i=0; i< len; i++) gsl_vector_set(v, i, 1./(i+1));
double square;
gsl_blas_ddot(v, v, &square);
#ifndef Testing
printf("1 + (1/2)^2 + (1/3)^2 + ...= %g\n", square);
#endif
double pi_over_six = gsl_pow_2(M_PI)/6.;
Diff(square, pi_over_six);
/* Now using apop_dot, in a few forms.
First, vector-as-data dot itself.
If one of the inputs is a vector,
apop_dot puts the output in a vector-as-data:*/
apop_data *v_as_data = apop_vector_to_data(v);
apop_data *vdotv = apop_dot(v_as_data, v_as_data);
Diff(gsl_vector_get(vdotv->vector, 0), pi_over_six);
/* As a matrix-as-data, via a compound literal.
Default is a len X 1 matrix. */
gsl_matrix *v_as_matrix = apop_vector_to_matrix(v);
apop_data dm = (apop_data){.matrix=v_as_matrix};
// (1 X len) dot (len X 1) --- OK, produce a scalar.
apop_data *mdotv = apop_dot(v_as_data, &dm);
double scalarval = gsl_vector_get(mdotv->vector, 0);
Diff(scalarval, pi_over_six);
//(len X 1) dot (len X 1) --- bad dimensions.
apop_opts.verbose=-1; //don't print an error.
apop_data *mdotv2 = apop_dot(&dm, v_as_data);
apop_opts.verbose=0; //back to sanity.
assert(mdotv2->error);
// If we want (len X 1) dot (1 X len) --> (len X len),
// use apop_vector_to_matrix.
apop_data dmr = (apop_data){.matrix=apop_vector_to_matrix(v, .row_col='r')};
apop_data *product_matrix = apop_dot(&dm, &dmr);
//The trace is the sum of squares:
gsl_vector_view trace = gsl_matrix_diagonal(product_matrix->matrix);
double tracesum = apop_sum(&trace.vector);
Diff(tracesum, pi_over_six);
apop_data_free(product_matrix);
gsl_matrix_free(dmr.matrix);
}
double apop_matrix_determinant ( const gsl_matrix *  in)

Find the determinant of a matrix. The in matrix is not destroyed in the process.

See also apop_matrix_inverse , or apop_det_and_inv to do both at once.

Parameters
inThe matrix to be determined.
Returns
The determinant.
gsl_matrix* apop_matrix_inverse ( const gsl_matrix *  in)

Inverts a matrix. The in matrix is not destroyed in the process. You may want to call apop_matrix_determinant first to check that your input is invertible, or use apop_det_and_inv to do both at once.

Parameters
inThe matrix to be inverted.
Returns
Its inverse.
apop_data* apop_matrix_pca ( gsl_matrix *  data,
int const  dimensions_we_want 
)

Principal component analysis: hand in a matrix and (optionally) a number of desired dimensions, and I'll return a data set where each column of the matrix is an eigenvector. The columns are sorted, so column zero has the greatest weight. The vector element of the data set gives the weights.

You also specify the number of elements your principal component space should have. If this is equal to the rank of the space in which the input data lives, then the sum of weights will be one. If the dimensions desired is less than that (probably so you can prepare a plot), then the weights will be accordingly smaller, giving you an indication of how much variation these dimensions explain.

Parameters
dataThe input matrix. (No default. If NULL, return NULL and print a warning iff apop_opts.verbose >= 1.) I modify int in place so that each column has mean zero.
dimensions_we_want(default: the size of the covariance matrix, i.e. data->size2) The singular value decomposition will return this many of the eigenvectors with the largest eigenvalues.
Returns
Returns a apop_data set whose matrix is the principal component space. Each column of the returned matrix will be another eigenvector; the columns will be ordered by the eigenvalues. The data set's vector will be the largest eigenvalues, scaled by the total of all eigenvalues (including those that were thrown out). The sum of these returned values will give you the percentage of variance explained by the factor analysis.
Exceptions
out->error=='a'Allocation error.
gsl_vector* apop_numerical_gradient ( apop_data data,
apop_model model,
double  delta 
)

The GSL provides one-dimensional numerical differentiation; here's the multidimensional extension.

Parameters
dataThe data set to use for all evaluations. It remains constant throughout.
modelThe model, expressing the function whose derivative is sought. The gradient is taken via small changes along the model parameters.
deltaThe size of the differential. If you explicitly give me a delta, I'll use it. If delta is not specified, but model has method_settings of type apop_ml_params, then the delta element is used for the differential. Else, I use 1e-3.
1 gsl_vector *gradient = apop_numerical_gradient(data, your_parametrized_model);

Autogenerated by doxygen on Sun Oct 26 2014 (Debian 0.999b+ds3-2).