function/************************************************************************* Weighted fitting by penalized cubic spline. Equidistant grid with M nodes on [min(x,xc),max(x,xc)] is used to build basis functions. Basis functions are cubic splines with natural boundary conditions. Problem is regularized by adding non-linearity penalty to the usual least squares penalty function: S(x) = arg min { LS + P }, where LS = SUM { w[i]^2*(y[i] - S(x[i]))^2 } - least squares penalty P = C*10^rho*integral{ S''(x)^2*dx } - non-linearity penalty rho - tunable constant given by user C - automatically determined scale parameter, makes penalty invariant with respect to scaling of X, Y, W. INPUT PARAMETERS: X - points, array[0..N-1]. Y - function values, array[0..N-1]. W - weights, array[0..N-1] Each summand in square sum of approximation deviations from given values is multiplied by the square of corresponding weight. Fill it by 1's if you don't want to solve weighted problem. N - number of points (optional): * N>0 * if given, only first N elements of X/Y/W are processed * if not given, automatically determined from X/Y/W sizes M - number of basis functions ( = number_of_nodes), M>=4. Rho - regularization constant passed by user. It penalizes nonlinearity in the regression spline. It is logarithmically scaled, i.e. actual value of regularization constant is calculated as 10^Rho. It is automatically scaled so that: * Rho=2.0 corresponds to moderate amount of nonlinearity * generally, it should be somewhere in the [-8.0,+8.0] If you do not want to penalize nonlineary, pass small Rho. Values as low as -15 should work. OUTPUT PARAMETERS: Info- same format as in LSFitLinearWC() subroutine. * Info>0 task is solved * Info<=0 an error occured: -4 means inconvergence of internal SVD or Cholesky decomposition; problem may be too ill-conditioned (very rare) S - spline interpolant. Rep - Following fields are set: * RMSError rms error on the (X,Y). * AvgError average error on the (X,Y). * AvgRelError average relative error on the non-zero Y * MaxError maximum error NON-WEIGHTED ERRORS ARE CALCULATED IMPORTANT: this subroitine doesn't calculate task's condition number for K<>0. NOTE 1: additional nodes are added to the spline outside of the fitting interval to force linearity when x<min(x,xc) or x>max(x,xc). It is done for consistency - we penalize non-linearity at [min(x,xc),max(x,xc)], so it is natural to force linearity outside of this interval. NOTE 2: function automatically sorts points, so caller may pass unsorted array. -- ALGLIB PROJECT -- Copyright 19.10.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.spline1dfitpenalizedw( double[] x, double[] y, double[] w, int m, double rho, out int info, out spline1dinterpolant s, out spline1dfitreport rep) public static void alglib.spline1dfitpenalizedw( double[] x, double[] y, double[] w, int n, int m, double rho, out int info, out spline1dinterpolant s, out spline1dfitreport rep)
Examples: [1]
public static int Main(string[] args) { // // In this example we demonstrate linear fitting by f(x|a) = a*exp(0.5*x). // // We have: // * y - vector of experimental data // * fmatrix - matrix of basis functions calculated at sample points // Actually, we have only one basis function F0 = exp(0.5*x). // double[,] fmatrix = new double[,]{{0.606531},{0.670320},{0.740818},{0.818731},{0.904837},{1.000000},{1.105171},{1.221403},{1.349859},{1.491825},{1.648721}}; double[] y = new double[]{1.133719,1.306522,1.504604,1.554663,1.884638,2.072436,2.257285,2.534068,2.622017,2.897713,3.219371}; int info; double[] c; alglib.lsfitreport rep; // // Linear fitting without weights // alglib.lsfitlinear(y, fmatrix, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(c,4)); // EXPECTED: [1.98650] // // Linear fitting with individual weights. // Slightly different result is returned. // double[] w = new double[]{1.414213,1,1,1,1,1,1,1,1,1,1}; alglib.lsfitlinearw(y, w, fmatrix, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(c,4)); // EXPECTED: [1.983354] System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // In this example we demonstrate linear fitting by f(x|a,b) = a*x+b // with simple constraint f(0)=0. // // We have: // * y - vector of experimental data // * fmatrix - matrix of basis functions sampled at [0,1] with step 0.2: // [ 1.0 0.0 ] // [ 1.0 0.2 ] // [ 1.0 0.4 ] // [ 1.0 0.6 ] // [ 1.0 0.8 ] // [ 1.0 1.0 ] // first column contains value of first basis function (constant term) // second column contains second basis function (linear term) // * cmatrix - matrix of linear constraints: // [ 1.0 0.0 0.0 ] // first two columns contain coefficients before basis functions, // last column contains desired value of their sum. // So [1,0,0] means "1*constant_term + 0*linear_term = 0" // double[] y = new double[]{0.072436,0.246944,0.491263,0.522300,0.714064,0.921929}; double[,] fmatrix = new double[,]{{1,0.0},{1,0.2},{1,0.4},{1,0.6},{1,0.8},{1,1.0}}; double[,] cmatrix = new double[,]{{1,0,0}}; int info; double[] c; alglib.lsfitreport rep; // // Constrained fitting without weights // alglib.lsfitlinearc(y, fmatrix, cmatrix, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(c,3)); // EXPECTED: [0,0.932933] // // Constrained fitting with individual weights // double[] w = new double[]{1,1.414213,1,1,1,1}; alglib.lsfitlinearwc(y, w, fmatrix, cmatrix, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(c,3)); // EXPECTED: [0,0.938322] System.Console.ReadLine(); return 0; }
public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)) // where x is a position on X-axis and c is adjustable parameter func = System.Math.Exp(-c[0]*x[0]*x[0]); } public static int Main(string[] args) { // // In this example we demonstrate exponential fitting // by f(x) = exp(-c*x^2) // using function value only. // // Gradient is estimated using combination of numerical differences // and secant updates. diffstep variable stores differentiation step // (we have to tell algorithm what step to use). // double[,] x = new double[,]{{-1},{-0.8},{-0.6},{-0.4},{-0.2},{0},{0.2},{0.4},{0.6},{0.8},{1.0}}; double[] y = new double[]{0.223130,0.382893,0.582748,0.786628,0.941765,1.000000,0.941765,0.786628,0.582748,0.382893,0.223130}; double[] c = new double[]{0.3}; double epsf = 0; double epsx = 0.000001; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; double diffstep = 0.0001; // // Fitting without weights // alglib.lsfitcreatef(x, y, c, diffstep, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.5] // // Fitting with weights // (you can change weights and see how it changes result) // double[] w = new double[]{1,1,1,1,1,1,1,1,1,1,1}; alglib.lsfitcreatewf(x, y, w, c, diffstep, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.5] System.Console.ReadLine(); return 0; }
public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)) // where x is a position on X-axis and c is adjustable parameter func = System.Math.Exp(-c[0]*x[0]*x[0]); } public static int Main(string[] args) { // // In this example we demonstrate exponential fitting by // f(x) = exp(-c*x^2) // subject to bound constraints // 0.0 <= c <= 1.0 // using function value only. // // Gradient is estimated using combination of numerical differences // and secant updates. diffstep variable stores differentiation step // (we have to tell algorithm what step to use). // // Unconstrained solution is c=1.5, but because of constraints we should // get c=1.0 (at the boundary). // double[,] x = new double[,]{{-1},{-0.8},{-0.6},{-0.4},{-0.2},{0},{0.2},{0.4},{0.6},{0.8},{1.0}}; double[] y = new double[]{0.223130,0.382893,0.582748,0.786628,0.941765,1.000000,0.941765,0.786628,0.582748,0.382893,0.223130}; double[] c = new double[]{0.3}; double[] bndl = new double[]{0.0}; double[] bndu = new double[]{1.0}; double epsf = 0; double epsx = 0.000001; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; double diffstep = 0.0001; alglib.lsfitcreatef(x, y, c, diffstep, out state); alglib.lsfitsetbc(state, bndl, bndu); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.0] System.Console.ReadLine(); return 0; }
public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)) // where x is a position on X-axis and c is adjustable parameter func = System.Math.Exp(-c[0]*x[0]*x[0]); } public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)) and gradient G={df/dc[i]} // where x is a position on X-axis and c is adjustable parameter. // IMPORTANT: gradient is calculated with respect to C, not to X func = System.Math.Exp(-c[0]*System.Math.Pow(x[0],2)); grad[0] = -System.Math.Pow(x[0],2)*func; } public static int Main(string[] args) { // // In this example we demonstrate exponential fitting // by f(x) = exp(-c*x^2) // using function value and gradient (with respect to c). // double[,] x = new double[,]{{-1},{-0.8},{-0.6},{-0.4},{-0.2},{0},{0.2},{0.4},{0.6},{0.8},{1.0}}; double[] y = new double[]{0.223130,0.382893,0.582748,0.786628,0.941765,1.000000,0.941765,0.786628,0.582748,0.382893,0.223130}; double[] c = new double[]{0.3}; double epsf = 0; double epsx = 0.000001; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; // // Fitting without weights // alglib.lsfitcreatefg(x, y, c, true, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.5] // // Fitting with weights // (you can change weights and see how it changes result) // double[] w = new double[]{1,1,1,1,1,1,1,1,1,1,1}; alglib.lsfitcreatewfg(x, y, w, c, true, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.5] System.Console.ReadLine(); return 0; }
public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)) // where x is a position on X-axis and c is adjustable parameter func = System.Math.Exp(-c[0]*x[0]*x[0]); } public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)) and gradient G={df/dc[i]} // where x is a position on X-axis and c is adjustable parameter. // IMPORTANT: gradient is calculated with respect to C, not to X func = System.Math.Exp(-c[0]*System.Math.Pow(x[0],2)); grad[0] = -System.Math.Pow(x[0],2)*func; } public static void function_cx_1_hess(double[] c, double[] x, ref double func, double[] grad, double[,] hess, object obj) { // this callback calculates f(c,x)=exp(-c0*sqr(x0)), gradient G={df/dc[i]} and Hessian H={d2f/(dc[i]*dc[j])} // where x is a position on X-axis and c is adjustable parameter. // IMPORTANT: gradient/Hessian are calculated with respect to C, not to X func = System.Math.Exp(-c[0]*System.Math.Pow(x[0],2)); grad[0] = -System.Math.Pow(x[0],2)*func; hess[0,0] = System.Math.Pow(x[0],4)*func; } public static int Main(string[] args) { // // In this example we demonstrate exponential fitting // by f(x) = exp(-c*x^2) // using function value, gradient and Hessian (with respect to c) // double[,] x = new double[,]{{-1},{-0.8},{-0.6},{-0.4},{-0.2},{0},{0.2},{0.4},{0.6},{0.8},{1.0}}; double[] y = new double[]{0.223130,0.382893,0.582748,0.786628,0.941765,1.000000,0.941765,0.786628,0.582748,0.382893,0.223130}; double[] c = new double[]{0.3}; double epsf = 0; double epsx = 0.000001; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; // // Fitting without weights // alglib.lsfitcreatefgh(x, y, c, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.5] // // Fitting with weights // (you can change weights and see how it changes result) // double[] w = new double[]{1,1,1,1,1,1,1,1,1,1,1}; alglib.lsfitcreatewfgh(x, y, w, c, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,1)); // EXPECTED: [1.5] System.Console.ReadLine(); return 0; }
public static void function_debt_func(double[] c, double[] x, ref double func, object obj) { // // this callback calculates f(c,x)=c[0]*(1+c[1]*(pow(x[0]-1999,c[2])-1)) // func = c[0]*(1+c[1]*(System.Math.Pow(x[0]-1999,c[2])-1)); } public static int Main(string[] args) { // // In this example we demonstrate fitting by // f(x) = c[0]*(1+c[1]*((x-1999)^c[2]-1)) // subject to bound constraints // -INF < c[0] < +INF // -10 <= c[1] <= +10 // 0.1 <= c[2] <= 2.0 // Data we want to fit are time series of Japan national debt // collected from 2000 to 2008 measured in USD (dollars, not // millions of dollars). // // Our variables are: // c[0] - debt value at initial moment (2000), // c[1] - direction coefficient (growth or decrease), // c[2] - curvature coefficient. // You may see that our variables are badly scaled - first one // is order of 10^12, and next two are somewhere about 1 in // magnitude. Such problem is difficult to solve without some // kind of scaling. // That is exactly where lsfitsetscale() function can be used. // We set scale of our variables to [1.0E12, 1, 1], which allows // us to easily solve this problem. // // You can try commenting out lsfitsetscale() call - and you will // see that algorithm will fail to converge. // double[,] x = new double[,]{{2000},{2001},{2002},{2003},{2004},{2005},{2006},{2007},{2008}}; double[] y = new double[]{4323239600000.0,4560913100000.0,5564091500000.0,6743189300000.0,7284064600000.0,7050129600000.0,7092221500000.0,8483907600000.0,8625804400000.0}; double[] c = new double[]{1.0e+13,1,1}; double epsf = 0; double epsx = 1.0e-5; double[] bndl = new double[]{-System.Double.PositiveInfinity,-10,0.1}; double[] bndu = new double[]{System.Double.PositiveInfinity,+10,2.0}; double[] s = new double[]{1.0e+12,1,1}; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; double diffstep = 1.0e-5; alglib.lsfitcreatef(x, y, c, diffstep, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitsetbc(state, bndl, bndu); alglib.lsfitsetscale(state, s); alglib.lsfitfit(state, function_debt_func, null, null); alglib.lsfitresults(state, out info, out c, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 2 System.Console.WriteLine("{0}", alglib.ap.format(c,-2)); // EXPECTED: [4.142560E+12, 0.434240, 0.565376] System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // This example demonstrates polynomial fitting. // // Fitting is done by two (M=2) functions from polynomial basis: // f0 = 1 // f1 = x // Basically, it just a linear fit; more complex polynomials may be used // (e.g. parabolas with M=3, cubic with M=4), but even such simple fit allows // us to demonstrate polynomialfit() function in action. // // We have: // * x set of abscissas // * y experimental data // // Additionally we demonstrate weighted fitting, where second point has // more weight than other ones. // double[] x = new double[]{0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}; double[] y = new double[]{0.00,0.05,0.26,0.32,0.33,0.43,0.60,0.60,0.77,0.98,1.02}; int m = 2; double t = 2; int info; alglib.barycentricinterpolant p; alglib.polynomialfitreport rep; double v; // // Fitting without individual weights // // NOTE: result is returned as barycentricinterpolant structure. // if you want to get representation in the power basis, // you can use barycentricbar2pow() function to convert // from barycentric to power representation (see docs for // POLINT subpackage for more info). // alglib.polynomialfit(x, y, m, out info, out p, out rep); v = alglib.barycentriccalc(p, t); System.Console.WriteLine("{0:F2", v); // EXPECTED: 2.011 // // Fitting with individual weights // // NOTE: slightly different result is returned // double[] w = new double[]{1,1.414213562,1,1,1,1,1,1,1,1,1}; double[] xc = new double[0]; double[] yc = new double[0]; int[] dc = new int[0]; alglib.polynomialfitwc(x, y, w, xc, yc, dc, m, out info, out p, out rep); v = alglib.barycentriccalc(p, t); System.Console.WriteLine("{0:F2", v); // EXPECTED: 2.023 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // This example demonstrates polynomial fitting. // // Fitting is done by two (M=2) functions from polynomial basis: // f0 = 1 // f1 = x // with simple constraint on function value // f(0) = 0 // Basically, it just a linear fit; more complex polynomials may be used // (e.g. parabolas with M=3, cubic with M=4), but even such simple fit allows // us to demonstrate polynomialfit() function in action. // // We have: // * x set of abscissas // * y experimental data // * xc points where constraints are placed // * yc constraints on derivatives // * dc derivative indices // (0 means function itself, 1 means first derivative) // double[] x = new double[]{1.0,1.0}; double[] y = new double[]{0.9,1.1}; double[] w = new double[]{1,1}; double[] xc = new double[]{0}; double[] yc = new double[]{0}; int[] dc = new int[]{0}; double t = 2; int m = 2; int info; alglib.barycentricinterpolant p; alglib.polynomialfitreport rep; double v; alglib.polynomialfitwc(x, y, w, xc, yc, dc, m, out info, out p, out rep); v = alglib.barycentriccalc(p, t); System.Console.WriteLine("{0:F2", v); // EXPECTED: 2.000 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // In this example we demonstrate penalized spline fitting of noisy data // // We have: // * x - abscissas // * y - vector of experimental data, straight line with small noise // double[] x = new double[]{0.00,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90}; double[] y = new double[]{0.10,0.00,0.30,0.40,0.30,0.40,0.62,0.68,0.75,0.95}; int info; double v; alglib.spline1dinterpolant s; alglib.spline1dfitreport rep; double rho; // // Fit with VERY small amount of smoothing (rho = -5.0) // and large number of basis functions (M=50).  // // With such regularization our spline should become close to the straight line fit. // We will compare its value in x=1.0 with results obtained from such fit. // rho = +10.0; alglib.spline1dfitpenalized(x, y, 50, rho, out info, out s, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 v = alglib.spline1dcalc(s, 1.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 0.969 // // In real life applications you may need some moderate degree of fitting, // so we try to fit once more with rho=3.0. // rho = +3.0; alglib.spline1dfitpenalized(x, y, 50, rho, out info, out s, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.ReadLine(); return 0; }
mannwhitneyu
subpackagemannwhitneyutest
function/************************************************************************* Mann-Whitney U-test This test checks hypotheses about whether X and Y are samples of two continuous distributions of the same shape and same median or whether their medians are different. The following tests are performed: * two-tailed test (null hypothesis - the medians are equal) * left-tailed test (null hypothesis - the median of the first sample is greater than or equal to the median of the second sample) * right-tailed test (null hypothesis - the median of the first sample is less than or equal to the median of the second sample). Requirements: * the samples are independent * X and Y are continuous distributions (or discrete distributions well- approximating continuous distributions) * distributions of X and Y have the same shape. The only possible difference is their position (i.e. the value of the median) * the number of elements in each sample is not less than 5 * the scale of measurement should be ordinal, interval or ratio (i.e. the test could not be applied to nominal variables). The test is non-parametric and doesn't require distributions to be normal. Input parameters: X - sample 1. Array whose index goes from 0 to N-1. N - size of the sample. N>=5 Y - sample 2. Array whose index goes from 0 to M-1. M - size of the sample. M>=5 Output parameters: BothTails - p-value for two-tailed test. If BothTails is less than the given significance level the null hypothesis is rejected. LeftTail - p-value for left-tailed test. If LeftTail is less than the given significance level, the null hypothesis is rejected. RightTail - p-value for right-tailed test. If RightTail is less than the given significance level the null hypothesis is rejected. To calculate p-values, special approximation is used. This method lets us calculate p-values with satisfactory accuracy in interval [0.0001, 1]. There is no approximation outside the [0.0001, 1] interval. Therefore, if the significance level outlies this interval, the test returns 0.0001. Relative precision of approximation of p-value: N M Max.err. Rms.err. 5..10 N..10 1.4e-02 6.0e-04 5..10 N..100 2.2e-02 5.3e-06 10..15 N..15 1.0e-02 3.2e-04 10..15 N..100 1.0e-02 2.2e-05 15..100 N..100 6.1e-03 2.7e-06 For N,M>100 accuracy checks weren't put into practice, but taking into account characteristics of asymptotic approximation used, precision should not be sharply different from the values for interval [5, 100]. -- ALGLIB -- Copyright 09.04.2007 by Bochkanov Sergey *************************************************************************/public static void alglib.mannwhitneyutest( double[] x, int n, double[] y, int m, out double bothtails, out double lefttail, out double righttail)
matdet
subpackage| matdet_d_1 | Determinant calculation, real matrix, short form | |
| matdet_d_2 | Determinant calculation, real matrix, full form | |
| matdet_d_3 | Determinant calculation, complex matrix, short form | |
| matdet_d_4 | Determinant calculation, complex matrix, full form | |
| matdet_d_5 | Determinant calculation, complex matrix with zero imaginary part, short form |
cmatrixdet
function/************************************************************************* Calculation of the determinant of a general matrix Input parameters: A - matrix, array[0..N-1, 0..N-1] N - (optional) size of matrix A: * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, automatically determined from matrix size (A must be square matrix) Result: determinant of matrix A. -- ALGLIB -- Copyright 2005 by Bochkanov Sergey *************************************************************************/public static complex alglib.cmatrixdet(complex[,] a) public static complex alglib.cmatrixdet(complex[,] a, int n)
cmatrixludet
function/************************************************************************* Determinant calculation of the matrix given by its LU decomposition. Input parameters: A - LU decomposition of the matrix (output of RMatrixLU subroutine). Pivots - table of permutations which were made during the LU decomposition. Output of RMatrixLU subroutine. N - (optional) size of matrix A: * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, automatically determined from matrix size (A must be square matrix) Result: matrix determinant. -- ALGLIB -- Copyright 2005 by Bochkanov Sergey *************************************************************************/public static complex alglib.cmatrixludet(complex[,] a, int[] pivots) public static complex alglib.cmatrixludet( complex[,] a, int[] pivots, int n)
rmatrixdet
function/************************************************************************* Calculation of the determinant of a general matrix Input parameters: A - matrix, array[0..N-1, 0..N-1] N - (optional) size of matrix A: * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, automatically determined from matrix size (A must be square matrix) Result: determinant of matrix A. -- ALGLIB -- Copyright 2005 by Bochkanov Sergey *************************************************************************/public static double alglib.rmatrixdet(double[,] a) public static double alglib.rmatrixdet(double[,] a, int n)
rmatrixludet
function/************************************************************************* Determinant calculation of the matrix given by its LU decomposition. Input parameters: A - LU decomposition of the matrix (output of RMatrixLU subroutine). Pivots - table of permutations which were made during the LU decomposition. Output of RMatrixLU subroutine. N - (optional) size of matrix A: * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, automatically determined from matrix size (A must be square matrix) Result: matrix determinant. -- ALGLIB -- Copyright 2005 by Bochkanov Sergey *************************************************************************/public static double alglib.rmatrixludet(double[,] a, int[] pivots) public static double alglib.rmatrixludet(double[,] a, int[] pivots, int n)
spdmatrixcholeskydet
function/************************************************************************* Determinant calculation of the matrix given by the Cholesky decomposition. Input parameters: A - Cholesky decomposition, output of SMatrixCholesky subroutine. N - (optional) size of matrix A: * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, automatically determined from matrix size (A must be square matrix) As the determinant is equal to the product of squares of diagonal elements, it’s not necessary to specify which triangle - lower or upper - the matrix is stored in. Result: matrix determinant. -- ALGLIB -- Copyright 2005-2008 by Bochkanov Sergey *************************************************************************/public static double alglib.spdmatrixcholeskydet(double[,] a) public static double alglib.spdmatrixcholeskydet(double[,] a, int n)
spdmatrixdet
function/************************************************************************* Determinant calculation of the symmetric positive definite matrix. Input parameters: A - matrix. Array with elements [0..N-1, 0..N-1]. N - (optional) size of matrix A: * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, automatically determined from matrix size (A must be square matrix) IsUpper - (optional) storage type: * if True, symmetric matrix A is given by its upper triangle, and the lower triangle isn’t used/changed by function * if False, symmetric matrix A is given by its lower triangle, and the upper triangle isn’t used/changed by function * if not given, both lower and upper triangles must be filled. Result: determinant of matrix A. If matrix A is not positive definite, exception is thrown. -- ALGLIB -- Copyright 2005-2008 by Bochkanov Sergey *************************************************************************/public static double alglib.spdmatrixdet(double[,] a) public static double alglib.spdmatrixdet(double[,] a, int n, bool isupper)
public static int Main(string[] args) { double[,] b = new double[,]{{1,2},{2,1}}; double a; a = alglib.rmatrixdet(b); System.Console.WriteLine("{0:F3", a); // EXPECTED: -3 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { double[,] b = new double[,]{{5,4},{4,5}}; double a; a = alglib.rmatrixdet(b, 2); System.Console.WriteLine("{0:F3", a); // EXPECTED: 9 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { alglib.complex[,] b = new alglib.complex[,]{{new alglib.complex(1,+1),2},{2,new alglib.complex(1,-1)}}; alglib.complex a; a = alglib.cmatrixdet(b); System.Console.WriteLine("{0}", alglib.ap.format(a,3)); // EXPECTED: -2 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { alglib.complex a; alglib.complex[,] b = new alglib.complex[,]{{new alglib.complex(0,5),4},{new alglib.complex(0,4),5}}; a = alglib.cmatrixdet(b, 2); System.Console.WriteLine("{0}", alglib.ap.format(a,3)); // EXPECTED: 9i System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { alglib.complex a; alglib.complex[,] b = new alglib.complex[,]{{9,1},{2,1}}; a = alglib.cmatrixdet(b); System.Console.WriteLine("{0}", alglib.ap.format(a,3)); // EXPECTED: 7 System.Console.ReadLine(); return 0; }
matgen
subpackagecmatrixrndcond
function/************************************************************************* Generation of random NxN complex matrix with given condition number C and norm2(A)=1 INPUT PARAMETERS: N - matrix size C - condition number (in 2-norm) OUTPUT PARAMETERS: A - random matrix with norm2(A)=1 and cond(A)=C -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixrndcond( int n, double c, out complex[,] a)
cmatrixrndorthogonal
function/************************************************************************* Generation of a random Haar distributed orthogonal complex matrix INPUT PARAMETERS: N - matrix size, N>=1 OUTPUT PARAMETERS: A - orthogonal NxN matrix, array[0..N-1,0..N-1] -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixrndorthogonal(int n, out complex[,] a)
cmatrixrndorthogonalfromtheleft
function/************************************************************************* Multiplication of MxN complex matrix by MxM random Haar distributed complex orthogonal matrix INPUT PARAMETERS: A - matrix, array[0..M-1, 0..N-1] M, N- matrix size OUTPUT PARAMETERS: A - Q*A, where Q is random MxM orthogonal matrix -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixrndorthogonalfromtheleft( ref complex[,] a, int m, int n)
cmatrixrndorthogonalfromtheright
function/************************************************************************* Multiplication of MxN complex matrix by NxN random Haar distributed complex orthogonal matrix INPUT PARAMETERS: A - matrix, array[0..M-1, 0..N-1] M, N- matrix size OUTPUT PARAMETERS: A - A*Q, where Q is random NxN orthogonal matrix -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixrndorthogonalfromtheright( ref complex[,] a, int m, int n)
hmatrixrndcond
function/************************************************************************* Generation of random NxN Hermitian matrix with given condition number and norm2(A)=1 INPUT PARAMETERS: N - matrix size C - condition number (in 2-norm) OUTPUT PARAMETERS: A - random matrix with norm2(A)=1 and cond(A)=C -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.hmatrixrndcond( int n, double c, out complex[,] a)
hmatrixrndmultiply
function/************************************************************************* Hermitian multiplication of NxN matrix by random Haar distributed complex orthogonal matrix INPUT PARAMETERS: A - matrix, array[0..N-1, 0..N-1] N - matrix size OUTPUT PARAMETERS: A - Q^H*A*Q, where Q is random NxN orthogonal matrix -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.hmatrixrndmultiply(ref complex[,] a, int n)
hpdmatrixrndcond
function/************************************************************************* Generation of random NxN Hermitian positive definite matrix with given condition number and norm2(A)=1 INPUT PARAMETERS: N - matrix size C - condition number (in 2-norm) OUTPUT PARAMETERS: A - random HPD matrix with norm2(A)=1 and cond(A)=C -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.hpdmatrixrndcond( int n, double c, out complex[,] a)
rmatrixrndcond
function/************************************************************************* Generation of random NxN matrix with given condition number and norm2(A)=1 INPUT PARAMETERS: N - matrix size C - condition number (in 2-norm) OUTPUT PARAMETERS: A - random matrix with norm2(A)=1 and cond(A)=C -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixrndcond(int n, double c, out double[,] a)
rmatrixrndorthogonal
function/************************************************************************* Generation of a random uniformly distributed (Haar) orthogonal matrix INPUT PARAMETERS: N - matrix size, N>=1 OUTPUT PARAMETERS: A - orthogonal NxN matrix, array[0..N-1,0..N-1] -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixrndorthogonal(int n, out double[,] a)
rmatrixrndorthogonalfromtheleft
function/************************************************************************* Multiplication of MxN matrix by MxM random Haar distributed orthogonal matrix INPUT PARAMETERS: A - matrix, array[0..M-1, 0..N-1] M, N- matrix size OUTPUT PARAMETERS: A - Q*A, where Q is random MxM orthogonal matrix -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixrndorthogonalfromtheleft( ref double[,] a, int m, int n)
rmatrixrndorthogonalfromtheright
function/************************************************************************* Multiplication of MxN matrix by NxN random Haar distributed orthogonal matrix INPUT PARAMETERS: A - matrix, array[0..M-1, 0..N-1] M, N- matrix size OUTPUT PARAMETERS: A - A*Q, where Q is random NxN orthogonal matrix -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixrndorthogonalfromtheright( ref double[,] a, int m, int n)
smatrixrndcond
function/************************************************************************* Generation of random NxN symmetric matrix with given condition number and norm2(A)=1 INPUT PARAMETERS: N - matrix size C - condition number (in 2-norm) OUTPUT PARAMETERS: A - random matrix with norm2(A)=1 and cond(A)=C -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.smatrixrndcond(int n, double c, out double[,] a)
smatrixrndmultiply
function/************************************************************************* Symmetric multiplication of NxN matrix by random Haar distributed orthogonal matrix INPUT PARAMETERS: A - matrix, array[0..N-1, 0..N-1] N - matrix size OUTPUT PARAMETERS: A - Q'*A*Q, where Q is random NxN orthogonal matrix -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.smatrixrndmultiply(ref double[,] a, int n)
spdmatrixrndcond
function/************************************************************************* Generation of random NxN symmetric positive definite matrix with given condition number and norm2(A)=1 INPUT PARAMETERS: N - matrix size C - condition number (in 2-norm) OUTPUT PARAMETERS: A - random SPD matrix with norm2(A)=1 and cond(A)=C -- ALGLIB routine -- 04.12.2009 Bochkanov Sergey *************************************************************************/public static void alglib.spdmatrixrndcond( int n, double c, out double[,] a)
matinv
subpackage| matinv_d_c1 | Complex matrix inverse | |
| matinv_d_hpd1 | HPD matrix inverse | |
| matinv_d_r1 | Real matrix inverse | |
| matinv_d_spd1 | SPD matrix inverse |
matinvreport
class/************************************************************************* Matrix inverse report: * R1 reciprocal of condition number in 1-norm * RInf reciprocal of condition number in inf-norm *************************************************************************/public class matinvreport { public double r1; public double rinf; }
cmatrixinverse
function/************************************************************************* Inversion of a general matrix. Input parameters: A - matrix N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) Output parameters: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse -- ALGLIB -- Copyright 2005 by Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixinverse( ref complex[,] a, out int info, out matinvreport rep) public static void alglib.cmatrixinverse( ref complex[,] a, int n, out int info, out matinvreport rep)
Examples: [1]
cmatrixluinverse
function/************************************************************************* Inversion of a matrix given by its LU decomposition. INPUT PARAMETERS: A - LU decomposition of the matrix (output of CMatrixLU subroutine). Pivots - table of permutations (the output of CMatrixLU subroutine). N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) OUTPUT PARAMETERS: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse -- ALGLIB routine -- 05.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixluinverse( ref complex[,] a, int[] pivots, out int info, out matinvreport rep) public static void alglib.cmatrixluinverse( ref complex[,] a, int[] pivots, int n, out int info, out matinvreport rep)
cmatrixtrinverse
function/************************************************************************* Triangular matrix inverse (complex) The subroutine inverts the following types of matrices: * upper triangular * upper triangular with unit diagonal * lower triangular * lower triangular with unit diagonal In case of an upper (lower) triangular matrix, the inverse matrix will also be upper (lower) triangular, and after the end of the algorithm, the inverse matrix replaces the source matrix. The elements below (above) the main diagonal are not changed by the algorithm. If the matrix has a unit diagonal, the inverse matrix also has a unit diagonal, and the diagonal elements are not passed to the algorithm. Input parameters: A - matrix, array[0..N-1, 0..N-1]. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) IsUpper - True, if the matrix is upper triangular. IsUnit - diagonal type (optional): * if True, matrix has unit diagonal (a[i,i] are NOT used) * if False, matrix diagonal is arbitrary * if not given, False is assumed Output parameters: Info - same as for RMatrixLUInverse Rep - same as for RMatrixLUInverse A - same as for RMatrixLUInverse. -- ALGLIB -- Copyright 05.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixtrinverse( ref complex[,] a, bool isupper, out int info, out matinvreport rep) public static void alglib.cmatrixtrinverse( ref complex[,] a, int n, bool isupper, bool isunit, out int info, out matinvreport rep)
hpdmatrixcholeskyinverse
function/************************************************************************* Inversion of a Hermitian positive definite matrix which is given by Cholesky decomposition. Input parameters: A - Cholesky decomposition of the matrix to be inverted: A=U’*U or A = L*L'. Output of HPDMatrixCholesky subroutine. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) IsUpper - storage type (optional): * if True, symmetric matrix A is given by its upper triangle, and the lower triangle isn’t used/changed by function * if False, symmetric matrix A is given by its lower triangle, and the upper triangle isn’t used/changed by function * if not given, lower half is used. Output parameters: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse -- ALGLIB routine -- 10.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.hpdmatrixcholeskyinverse( ref complex[,] a, out int info, out matinvreport rep) public static void alglib.hpdmatrixcholeskyinverse( ref complex[,] a, int n, bool isupper, out int info, out matinvreport rep)
hpdmatrixinverse
function/************************************************************************* Inversion of a Hermitian positive definite matrix. Given an upper or lower triangle of a Hermitian positive definite matrix, the algorithm generates matrix A^-1 and saves the upper or lower triangle depending on the input. Input parameters: A - matrix to be inverted (upper or lower triangle). Array with elements [0..N-1,0..N-1]. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) IsUpper - storage type (optional): * if True, symmetric matrix A is given by its upper triangle, and the lower triangle isn’t used/changed by function * if False, symmetric matrix A is given by its lower triangle, and the upper triangle isn’t used/changed by function * if not given, both lower and upper triangles must be filled. Output parameters: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse -- ALGLIB routine -- 10.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.hpdmatrixinverse( ref complex[,] a, out int info, out matinvreport rep) public static void alglib.hpdmatrixinverse( ref complex[,] a, int n, bool isupper, out int info, out matinvreport rep)
Examples: [1]
rmatrixinverse
function/************************************************************************* Inversion of a general matrix. Input parameters: A - matrix. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) Output parameters: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse Result: True, if the matrix is not singular. False, if the matrix is singular. -- ALGLIB -- Copyright 2005-2010 by Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixinverse( ref double[,] a, out int info, out matinvreport rep) public static void alglib.rmatrixinverse( ref double[,] a, int n, out int info, out matinvreport rep)
Examples: [1]
rmatrixluinverse
function/************************************************************************* Inversion of a matrix given by its LU decomposition. INPUT PARAMETERS: A - LU decomposition of the matrix (output of RMatrixLU subroutine). Pivots - table of permutations (the output of RMatrixLU subroutine). N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) OUTPUT PARAMETERS: Info - return code: * -3 A is singular, or VERY close to singular. it is filled by zeros in such cases. * 1 task is solved (but matrix A may be ill-conditioned, check R1/RInf parameters for condition numbers). Rep - solver report, see below for more info A - inverse of matrix A. Array whose indexes range within [0..N-1, 0..N-1]. SOLVER REPORT Subroutine sets following fields of the Rep structure: * R1 reciprocal of condition number: 1/cond(A), 1-norm. * RInf reciprocal of condition number: 1/cond(A), inf-norm. -- ALGLIB routine -- 05.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixluinverse( ref double[,] a, int[] pivots, out int info, out matinvreport rep) public static void alglib.rmatrixluinverse( ref double[,] a, int[] pivots, int n, out int info, out matinvreport rep)
rmatrixtrinverse
function/************************************************************************* Triangular matrix inverse (real) The subroutine inverts the following types of matrices: * upper triangular * upper triangular with unit diagonal * lower triangular * lower triangular with unit diagonal In case of an upper (lower) triangular matrix, the inverse matrix will also be upper (lower) triangular, and after the end of the algorithm, the inverse matrix replaces the source matrix. The elements below (above) the main diagonal are not changed by the algorithm. If the matrix has a unit diagonal, the inverse matrix also has a unit diagonal, and the diagonal elements are not passed to the algorithm. Input parameters: A - matrix, array[0..N-1, 0..N-1]. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) IsUpper - True, if the matrix is upper triangular. IsUnit - diagonal type (optional): * if True, matrix has unit diagonal (a[i,i] are NOT used) * if False, matrix diagonal is arbitrary * if not given, False is assumed Output parameters: Info - same as for RMatrixLUInverse Rep - same as for RMatrixLUInverse A - same as for RMatrixLUInverse. -- ALGLIB -- Copyright 05.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixtrinverse( ref double[,] a, bool isupper, out int info, out matinvreport rep) public static void alglib.rmatrixtrinverse( ref double[,] a, int n, bool isupper, bool isunit, out int info, out matinvreport rep)
spdmatrixcholeskyinverse
function/************************************************************************* Inversion of a symmetric positive definite matrix which is given by Cholesky decomposition. Input parameters: A - Cholesky decomposition of the matrix to be inverted: A=U’*U or A = L*L'. Output of SPDMatrixCholesky subroutine. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) IsUpper - storage type (optional): * if True, symmetric matrix A is given by its upper triangle, and the lower triangle isn’t used/changed by function * if False, symmetric matrix A is given by its lower triangle, and the upper triangle isn’t used/changed by function * if not given, lower half is used. Output parameters: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse -- ALGLIB routine -- 10.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.spdmatrixcholeskyinverse( ref double[,] a, out int info, out matinvreport rep) public static void alglib.spdmatrixcholeskyinverse( ref double[,] a, int n, bool isupper, out int info, out matinvreport rep)
spdmatrixinverse
function/************************************************************************* Inversion of a symmetric positive definite matrix. Given an upper or lower triangle of a symmetric positive definite matrix, the algorithm generates matrix A^-1 and saves the upper or lower triangle depending on the input. Input parameters: A - matrix to be inverted (upper or lower triangle). Array with elements [0..N-1,0..N-1]. N - size of matrix A (optional) : * if given, only principal NxN submatrix is processed and overwritten. other elements are unchanged. * if not given, size is automatically determined from matrix size (A must be square matrix) IsUpper - storage type (optional): * if True, symmetric matrix A is given by its upper triangle, and the lower triangle isn’t used/changed by function * if False, symmetric matrix A is given by its lower triangle, and the upper triangle isn’t used/changed by function * if not given, both lower and upper triangles must be filled. Output parameters: Info - return code, same as in RMatrixLUInverse Rep - solver report, same as in RMatrixLUInverse A - inverse of matrix A, same as in RMatrixLUInverse -- ALGLIB routine -- 10.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.spdmatrixinverse( ref double[,] a, out int info, out matinvreport rep) public static void alglib.spdmatrixinverse( ref double[,] a, int n, bool isupper, out int info, out matinvreport rep)
Examples: [1]
public static int Main(string[] args) { alglib.complex[,] a = new alglib.complex[,]{{new alglib.complex(0,1),-1},{new alglib.complex(0,1),1}}; int info; alglib.matinvreport rep; alglib.cmatrixinverse(ref a, out info, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(a,4)); // EXPECTED: [[-0.5i,-0.5i],[-0.5,0.5]] System.Console.WriteLine("{0:F4", rep.r1); // EXPECTED: 0.5 System.Console.WriteLine("{0:F4", rep.rinf); // EXPECTED: 0.5 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { alglib.complex[,] a = new alglib.complex[,]{{2,1},{1,2}}; int info; alglib.matinvreport rep; alglib.hpdmatrixinverse(ref a, out info, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(a,4)); // EXPECTED: [[0.666666,-0.333333],[-0.333333,0.666666]] System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { double[,] a = new double[,]{{1,-1},{1,1}}; int info; alglib.matinvreport rep; alglib.rmatrixinverse(ref a, out info, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(a,4)); // EXPECTED: [[0.5,0.5],[-0.5,0.5]] System.Console.WriteLine("{0:F4", rep.r1); // EXPECTED: 0.5 System.Console.WriteLine("{0:F4", rep.rinf); // EXPECTED: 0.5 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { double[,] a = new double[,]{{2,1},{1,2}}; int info; alglib.matinvreport rep; alglib.spdmatrixinverse(ref a, out info, out rep); System.Console.WriteLine("{0}", info); // EXPECTED: 1 System.Console.WriteLine("{0}", alglib.ap.format(a,4)); // EXPECTED: [[0.666666,-0.333333],[-0.333333,0.666666]] System.Console.ReadLine(); return 0; }
mcpd
subpackage| mcpd_simple1 | Simple unconstrained MCPD model (no entry/exit states) | |
| mcpd_simple2 | Simple MCPD model (no entry/exit states) with equality constraints |
mcpdreport
class/************************************************************************* This structure is a MCPD training report: InnerIterationsCount - number of inner iterations of the underlying optimization algorithm OuterIterationsCount - number of outer iterations of the underlying optimization algorithm NFEV - number of merit function evaluations TerminationType - termination type (same as for MinBLEIC optimizer, positive values denote success, negative ones - failure) -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public class mcpdreport { public int inneriterationscount; public int outeriterationscount; public int nfev; public int terminationtype; }
mcpdstate
class/************************************************************************* This structure is a MCPD (Markov Chains for Population Data) solver. You should use ALGLIB functions in order to work with this object. -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public class mcpdstate { }
mcpdaddbc
function/************************************************************************* This function is used to add bound constraints on the elements of the transition matrix P. MCPD solver has four types of constraints which can be placed on P: * user-specified equality constraints (optional) * user-specified bound constraints (optional) * user-specified general linear constraints (optional) * basic constraints (always present): * non-negativity: P[i,j]>=0 * consistency: every column of P sums to 1.0 Final constraints which are passed to the underlying optimizer are calculated as intersection of all present constraints. For example, you may specify boundary constraint on P[0,0] and equality one: 0.1<=P[0,0]<=0.9 P[0,0]=0.5 Such combination of constraints will be silently reduced to their intersection, which is P[0,0]=0.5. This function can be used to ADD bound constraint for one element of P without changing constraints for other elements. You can also use MCPDSetBC() function which allows to place bound constraints on arbitrary subset of elements of P. Set of constraints is specified by BndL/BndU matrices, which may contain arbitrary combination of finite numbers or infinities (like -INF<x<=0.5 or 0.1<=x<+INF). These functions (MCPDSetBC and MCPDAddBC) interact as follows: * there is internal matrix of bound constraints which is stored in the MCPD solver * MCPDSetBC() replaces this matrix by another one (SET) * MCPDAddBC() modifies one element of this matrix and leaves other ones unchanged (ADD) * thus MCPDAddBC() call preserves all modifications done by previous calls, while MCPDSetBC() completely discards all changes done to the equality constraints. INPUT PARAMETERS: S - solver I - row index of element being constrained J - column index of element being constrained BndL - lower bound BndU - upper bound -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdaddbc( mcpdstate s, int i, int j, double bndl, double bndu)
mcpdaddec
function/************************************************************************* This function is used to add equality constraints on the elements of the transition matrix P. MCPD solver has four types of constraints which can be placed on P: * user-specified equality constraints (optional) * user-specified bound constraints (optional) * user-specified general linear constraints (optional) * basic constraints (always present): * non-negativity: P[i,j]>=0 * consistency: every column of P sums to 1.0 Final constraints which are passed to the underlying optimizer are calculated as intersection of all present constraints. For example, you may specify boundary constraint on P[0,0] and equality one: 0.1<=P[0,0]<=0.9 P[0,0]=0.5 Such combination of constraints will be silently reduced to their intersection, which is P[0,0]=0.5. This function can be used to ADD equality constraint for one element of P without changing constraints for other elements. You can also use MCPDSetEC() function which allows you to specify arbitrary set of equality constraints in one call. These functions (MCPDSetEC and MCPDAddEC) interact as follows: * there is internal matrix of equality constraints which is stored in the MCPD solver * MCPDSetEC() replaces this matrix by another one (SET) * MCPDAddEC() modifies one element of this matrix and leaves other ones unchanged (ADD) * thus MCPDAddEC() call preserves all modifications done by previous calls, while MCPDSetEC() completely discards all changes done to the equality constraints. INPUT PARAMETERS: S - solver I - row index of element being constrained J - column index of element being constrained C - value (constraint for P[I,J]). Can be either NAN (no constraint) or finite value from [0,1]. NOTES: 1. infinite values of C will lead to exception being thrown. Values less than 0.0 or greater than 1.0 will lead to error code being returned after call to MCPDSolve(). -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdaddec(mcpdstate s, int i, int j, double c)
Examples: [1]
mcpdaddtrack
function/************************************************************************* This function is used to add a track - sequence of system states at the different moments of its evolution. You may add one or several tracks to the MCPD solver.  For example, if you pass two tracks, A1-A2-A3 (system at t=A+1, t=A+2 and t=A+3) and B1-B2-B3, then solver will try to model transitions from t=A+1 to t=A+2, t=A+2 to t=A+3, t=B+1 to t=B+2, t=B+2 to t=B+3. But it WONT mix these two tracks - i.e. it wont try to model transition from t=A+3 to t=B+1. INPUT PARAMETERS: S - solver XY - track, array[K,N]: * I-th row is a state at t=I * elements of XY must be non-negative (exception will be thrown on negative elements) K - number of points in a track * if given, only leading K rows of XY are used * if not given, automatically determined from size of XY NOTES: 1. Track may contain either proportional or population data: * with proportional data all rows of XY must sum to 1.0, i.e. we have proportions instead of absolute population values * with population data rows of XY contain population counts and generally do not sum to 1.0 (although they still must be non-negative) -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdaddtrack(mcpdstate s, double[,] xy) public static void alglib.mcpdaddtrack(mcpdstate s, double[,] xy, int k)
mcpdcreate
function/************************************************************************* DESCRIPTION: This function creates MCPD (Markov Chains for Population Data) solver. This solver can be used to find transition matrix P for N-dimensional prediction problem where transition from X[i] to X[i+1] is modelled as X[i+1] = P*X[i] where X[i] and X[i+1] are N-dimensional population vectors (components of each X are non-negative), and P is a N*N transition matrix (elements of P are non-negative, each column sums to 1.0). Such models arise when when: * there is some population of individuals * individuals can have different states * individuals can transit from one state to another * population size is constant, i.e. there is no new individuals and no one leaves population * you want to model transitions of individuals from one state into another USAGE: Here we give very brief outline of the MCPD. We strongly recommend you to read examples in the ALGLIB Reference Manual and to read ALGLIB User Guide on data analysis which is available at http://www.alglib.net/dataanalysis/ 1. User initializes algorithm state with MCPDCreate() call 2. User adds one or more tracks - sequences of states which describe evolution of a system being modelled from different starting conditions 3. User may add optional boundary, equality and/or linear constraints on the coefficients of P by calling one of the following functions: * MCPDSetEC() to set equality constraints * MCPDSetBC() to set bound constraints * MCPDSetLC() to set linear constraints 4. Optionally, user may set custom weights for prediction errors (by default, algorithm assigns non-equal, automatically chosen weights for errors in the prediction of different components of X). It can be done with a call of MCPDSetPredictionWeights() function. 5. User calls MCPDSolve() function which takes algorithm state and pointer (delegate, etc.) to callback function which calculates F/G. 6. User calls MCPDResults() to get solution INPUT PARAMETERS: N - problem dimension, N>=1 OUTPUT PARAMETERS: State - structure stores algorithm state -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdcreate(int n, out mcpdstate s)
mcpdcreateentry
function/************************************************************************* DESCRIPTION: This function is a specialized version of MCPDCreate() function, and we recommend you to read comments for this function for general information about MCPD solver. This function creates MCPD (Markov Chains for Population Data) solver for "Entry-state" model, i.e. model where transition from X[i] to X[i+1] is modelled as X[i+1] = P*X[i] where X[i] and X[i+1] are N-dimensional state vectors P is a N*N transition matrix and one selected component of X[] is called "entry" state and is treated in a special way: system state always transits from "entry" state to some another state system state can not transit from any state into "entry" state Such conditions basically mean that row of P which corresponds to "entry" state is zero. Such models arise when: * there is some population of individuals * individuals can have different states * individuals can transit from one state to another * population size is NOT constant - at every moment of time there is some (unpredictable) amount of "new" individuals, which can transit into one of the states at the next turn, but still no one leaves population * you want to model transitions of individuals from one state into another * but you do NOT want to predict amount of "new" individuals because it does not depends on individuals already present (hence system can not transit INTO entry state - it can only transit FROM it). This model is discussed in more details in the ALGLIB User Guide (see http://www.alglib.net/dataanalysis/ for more data). INPUT PARAMETERS: N - problem dimension, N>=2 EntryState- index of entry state, in 0..N-1 OUTPUT PARAMETERS: State - structure stores algorithm state -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdcreateentry( int n, int entrystate, out mcpdstate s)
mcpdcreateentryexit
function/************************************************************************* DESCRIPTION: This function is a specialized version of MCPDCreate() function, and we recommend you to read comments for this function for general information about MCPD solver. This function creates MCPD (Markov Chains for Population Data) solver for "Entry-Exit-states" model, i.e. model where transition from X[i] to X[i+1] is modelled as X[i+1] = P*X[i] where X[i] and X[i+1] are N-dimensional state vectors P is a N*N transition matrix one selected component of X[] is called "entry" state and is treated in a special way: system state always transits from "entry" state to some another state system state can not transit from any state into "entry" state and another one component of X[] is called "exit" state and is treated in a special way too: system state can transit from any state into "exit" state system state can not transit from "exit" state into any other state transition operator discards "exit" state (makes it zero at each turn) Such conditions basically mean that: row of P which corresponds to "entry" state is zero column of P which corresponds to "exit" state is zero Multiplication by such P may decrease sum of vector components. Such models arise when: * there is some population of individuals * individuals can have different states * individuals can transit from one state to another * population size is NOT constant * at every moment of time there is some (unpredictable) amount of "new" individuals, which can transit into one of the states at the next turn * some individuals can move (predictably) into "exit" state and leave population at the next turn * you want to model transitions of individuals from one state into another, including transitions from the "entry" state and into the "exit" state. * but you do NOT want to predict amount of "new" individuals because it does not depends on individuals already present (hence system can not transit INTO entry state - it can only transit FROM it). This model is discussed in more details in the ALGLIB User Guide (see http://www.alglib.net/dataanalysis/ for more data). INPUT PARAMETERS: N - problem dimension, N>=2 EntryState- index of entry state, in 0..N-1 ExitState- index of exit state, in 0..N-1 OUTPUT PARAMETERS: State - structure stores algorithm state -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdcreateentryexit( int n, int entrystate, int exitstate, out mcpdstate s)
mcpdcreateexit
function/************************************************************************* DESCRIPTION: This function is a specialized version of MCPDCreate() function, and we recommend you to read comments for this function for general information about MCPD solver. This function creates MCPD (Markov Chains for Population Data) solver for "Exit-state" model, i.e. model where transition from X[i] to X[i+1] is modelled as X[i+1] = P*X[i] where X[i] and X[i+1] are N-dimensional state vectors P is a N*N transition matrix and one selected component of X[] is called "exit" state and is treated in a special way: system state can transit from any state into "exit" state system state can not transit from "exit" state into any other state transition operator discards "exit" state (makes it zero at each turn) Such conditions basically mean that column of P which corresponds to "exit" state is zero. Multiplication by such P may decrease sum of vector components. Such models arise when: * there is some population of individuals * individuals can have different states * individuals can transit from one state to another * population size is NOT constant - individuals can move into "exit" state and leave population at the next turn, but there are no new individuals * amount of individuals which leave population can be predicted * you want to model transitions of individuals from one state into another (including transitions into the "exit" state) This model is discussed in more details in the ALGLIB User Guide (see http://www.alglib.net/dataanalysis/ for more data). INPUT PARAMETERS: N - problem dimension, N>=2 ExitState- index of exit state, in 0..N-1 OUTPUT PARAMETERS: State - structure stores algorithm state -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdcreateexit( int n, int exitstate, out mcpdstate s)
mcpdresults
function/************************************************************************* MCPD results INPUT PARAMETERS: State - algorithm state OUTPUT PARAMETERS: P - array[N,N], transition matrix Rep - optimization report. You should check Rep.TerminationType in order to distinguish successful termination from unsuccessful one. Speaking short, positive values denote success, negative ones are failures. More information about fields of this structure can be found in the comments on MCPDReport datatype. -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdresults( mcpdstate s, out double[,] p, out mcpdreport rep)
mcpdsetbc
function/************************************************************************* This function is used to add bound constraints on the elements of the transition matrix P. MCPD solver has four types of constraints which can be placed on P: * user-specified equality constraints (optional) * user-specified bound constraints (optional) * user-specified general linear constraints (optional) * basic constraints (always present): * non-negativity: P[i,j]>=0 * consistency: every column of P sums to 1.0 Final constraints which are passed to the underlying optimizer are calculated as intersection of all present constraints. For example, you may specify boundary constraint on P[0,0] and equality one: 0.1<=P[0,0]<=0.9 P[0,0]=0.5 Such combination of constraints will be silently reduced to their intersection, which is P[0,0]=0.5. This function can be used to place bound constraints on arbitrary subset of elements of P. Set of constraints is specified by BndL/BndU matrices, which may contain arbitrary combination of finite numbers or infinities (like -INF<x<=0.5 or 0.1<=x<+INF). You can also use MCPDAddBC() function which allows to ADD bound constraint for one element of P without changing constraints for other elements. These functions (MCPDSetBC and MCPDAddBC) interact as follows: * there is internal matrix of bound constraints which is stored in the MCPD solver * MCPDSetBC() replaces this matrix by another one (SET) * MCPDAddBC() modifies one element of this matrix and leaves other ones unchanged (ADD) * thus MCPDAddBC() call preserves all modifications done by previous calls, while MCPDSetBC() completely discards all changes done to the equality constraints. INPUT PARAMETERS: S - solver BndL - lower bounds constraints, array[N,N]. Elements of BndL can be finite numbers or -INF. BndU - upper bounds constraints, array[N,N]. Elements of BndU can be finite numbers or +INF. -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsetbc( mcpdstate s, double[,] bndl, double[,] bndu)
mcpdsetec
function/************************************************************************* This function is used to add equality constraints on the elements of the transition matrix P. MCPD solver has four types of constraints which can be placed on P: * user-specified equality constraints (optional) * user-specified bound constraints (optional) * user-specified general linear constraints (optional) * basic constraints (always present): * non-negativity: P[i,j]>=0 * consistency: every column of P sums to 1.0 Final constraints which are passed to the underlying optimizer are calculated as intersection of all present constraints. For example, you may specify boundary constraint on P[0,0] and equality one: 0.1<=P[0,0]<=0.9 P[0,0]=0.5 Such combination of constraints will be silently reduced to their intersection, which is P[0,0]=0.5. This function can be used to place equality constraints on arbitrary subset of elements of P. Set of constraints is specified by EC, which may contain either NAN's or finite numbers from [0,1]. NAN denotes absence of constraint, finite number denotes equality constraint on specific element of P. You can also use MCPDAddEC() function which allows to ADD equality constraint for one element of P without changing constraints for other elements. These functions (MCPDSetEC and MCPDAddEC) interact as follows: * there is internal matrix of equality constraints which is stored in the MCPD solver * MCPDSetEC() replaces this matrix by another one (SET) * MCPDAddEC() modifies one element of this matrix and leaves other ones unchanged (ADD) * thus MCPDAddEC() call preserves all modifications done by previous calls, while MCPDSetEC() completely discards all changes done to the equality constraints. INPUT PARAMETERS: S - solver EC - equality constraints, array[N,N]. Elements of EC can be either NAN's or finite numbers from [0,1]. NAN denotes absence of constraints, while finite value denotes equality constraint on the corresponding element of P. NOTES: 1. infinite values of EC will lead to exception being thrown. Values less than 0.0 or greater than 1.0 will lead to error code being returned after call to MCPDSolve(). -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsetec(mcpdstate s, double[,] ec)
mcpdsetlc
function/************************************************************************* This function is used to set linear equality/inequality constraints on the elements of the transition matrix P. This function can be used to set one or several general linear constraints on the elements of P. Two types of constraints are supported: * equality constraints * inequality constraints (both less-or-equal and greater-or-equal) Coefficients of constraints are specified by matrix C (one of the parameters). One row of C corresponds to one constraint. Because transition matrix P has N*N elements, we need N*N columns to store all coefficients (they are stored row by row), and one more column to store right part - hence C has N*N+1 columns. Constraint kind is stored in the CT array. Thus, I-th linear constraint is P[0,0]*C[I,0] + P[0,1]*C[I,1] + .. + P[0,N-1]*C[I,N-1] + + P[1,0]*C[I,N] + P[1,1]*C[I,N+1] + ... + + P[N-1,N-1]*C[I,N*N-1] ?=? C[I,N*N] where ?=? can be either "=" (CT[i]=0), "<=" (CT[i]<0) or ">=" (CT[i]>0). Your constraint may involve only some subset of P (less than N*N elements). For example it can be something like P[0,0] + P[0,1] = 0.5 In this case you still should pass matrix with N*N+1 columns, but all its elements (except for C[0,0], C[0,1] and C[0,N*N-1]) will be zero. INPUT PARAMETERS: S - solver C - array[K,N*N+1] - coefficients of constraints (see above for complete description) CT - array[K] - constraint types (see above for complete description) K - number of equality/inequality constraints, K>=0: * if given, only leading K elements of C/CT are used * if not given, automatically determined from sizes of C/CT -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsetlc(mcpdstate s, double[,] c, int[] ct) public static void alglib.mcpdsetlc( mcpdstate s, double[,] c, int[] ct, int k)
mcpdsetpredictionweights
function/************************************************************************* This function is used to change prediction weights MCPD solver scales prediction errors as follows Error(P) = ||W*(y-P*x)||^2 where x is a system state at time t y is a system state at time t+1 P is a transition matrix W is a diagonal scaling matrix By default, weights are chosen in order to minimize relative prediction error instead of absolute one. For example, if one component of state is about 0.5 in magnitude and another one is about 0.05, then algorithm will make corresponding weights equal to 2.0 and 20.0. INPUT PARAMETERS: S - solver PW - array[N], weights: * must be non-negative values (exception will be thrown otherwise) * zero values will be replaced by automatically chosen values -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsetpredictionweights( mcpdstate s, double[] pw)
mcpdsetprior
function/************************************************************************* This function allows to set prior values used for regularization of your problem. By default, regularizing term is equal to r*||P-prior_P||^2, where r is a small non-zero value, P is transition matrix, prior_P is identity matrix, ||X||^2 is a sum of squared elements of X. This function allows you to change prior values prior_P. You can also change r with MCPDSetTikhonovRegularizer() function. INPUT PARAMETERS: S - solver PP - array[N,N], matrix of prior values: 1. elements must be real numbers from [0,1] 2. columns must sum to 1.0. First property is checked (exception is thrown otherwise), while second one is not checked/enforced. -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsetprior(mcpdstate s, double[,] pp)
mcpdsettikhonovregularizer
function/************************************************************************* This function allows to tune amount of Tikhonov regularization being applied to your problem. By default, regularizing term is equal to r*||P-prior_P||^2, where r is a small non-zero value, P is transition matrix, prior_P is identity matrix, ||X||^2 is a sum of squared elements of X. This function allows you to change coefficient r. You can also change prior values with MCPDSetPrior() function. INPUT PARAMETERS: S - solver V - regularization coefficient, finite non-negative value. It is not recommended to specify zero value unless you are pretty sure that you want it. -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsettikhonovregularizer( mcpdstate s, double v)
mcpdsolve
function/************************************************************************* This function is used to start solution of the MCPD problem. After return from this function, you can use MCPDResults() to get solution and completion code. -- ALGLIB -- Copyright 23.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mcpdsolve(mcpdstate s)
public static int Main(string[] args) { // // The very simple MCPD example // // We have a loan portfolio. Our loans can be in one of two states: // * normal loans ("good" ones) // * past due loans ("bad" ones) // // We assume that: // * loans can transition from any state to any other state. In // particular, past due loan can become "good" one at any moment // with same (fixed) probability. Not realistic, but it is toy example :) // * portfolio size does not change over time // // Thus, we have following model // state_new = P*state_old // where // ( p00 p01 ) // P = ( ) // ( p10 p11 ) // // We want to model transitions between these two states using MCPD // approach (Markov Chains for Proportional/Population Data), i.e. // to restore hidden transition matrix P using actual portfolio data. // We have: // * poportional data, i.e. proportion of loans in the normal and past // due states (not portfolio size measured in some currency, although // it is possible to work with population data too) // * two tracks, i.e. two sequences which describe portfolio // evolution from two different starting states: [1,0] (all loans // are "good") and [0.8,0.2] (only 80% of portfolio is in the "good" // state) // alglib.mcpdstate s; alglib.mcpdreport rep; double[,] p; double[,] track0 = new double[,]{{1.00000,0.00000},{0.95000,0.05000},{0.92750,0.07250},{0.91738,0.08263},{0.91282,0.08718}}; double[,] track1 = new double[,]{{0.80000,0.20000},{0.86000,0.14000},{0.88700,0.11300},{0.89915,0.10085}}; alglib.mcpdcreate(2, out s); alglib.mcpdaddtrack(s, track0); alglib.mcpdaddtrack(s, track1); alglib.mcpdsolve(s); alglib.mcpdresults(s, out p, out rep); // // Hidden matrix P is equal to // ( 0.95 0.50 ) // ( ) // ( 0.05 0.50 ) // which means that "good" loans can become "bad" with 5% probability, // while "bad" loans will return to good state with 50% probability. // System.Console.WriteLine("{0}", alglib.ap.format(p,2)); // EXPECTED: [[0.95,0.50],[0.05,0.50]] System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // Simple MCPD example // // We have a loan portfolio. Our loans can be in one of three states: // * normal loans // * past due loans // * charged off loans // // We assume that: // * normal loan can stay normal or become past due (but not charged off) // * past due loan can stay past due, become normal or charged off // * charged off loan will stay charged off for the rest of eternity // * portfolio size does not change over time // Not realistic, but it is toy example :) // // Thus, we have following model // state_new = P*state_old // where // ( p00 p01 ) // P = ( p10 p11 ) // ( p21 1 ) // i.e. four elements of P are known a priori. // // Although it is possible (given enough data) to In order to enforce // this property we set equality constraints on these elements. // // We want to model transitions between these two states using MCPD // approach (Markov Chains for Proportional/Population Data), i.e. // to restore hidden transition matrix P using actual portfolio data. // We have: // * poportional data, i.e. proportion of loans in the current and past // due states (not portfolio size measured in some currency, although // it is possible to work with population data too) // * two tracks, i.e. two sequences which describe portfolio // evolution from two different starting states: [1,0,0] (all loans // are "good") and [0.8,0.2,0.0] (only 80% of portfolio is in the "good" // state) // alglib.mcpdstate s; alglib.mcpdreport rep; double[,] p; double[,] track0 = new double[,]{{1.000000,0.000000,0.000000},{0.950000,0.050000,0.000000},{0.927500,0.060000,0.012500},{0.911125,0.061375,0.027500},{0.896256,0.060900,0.042844}}; double[,] track1 = new double[,]{{0.800000,0.200000,0.000000},{0.860000,0.090000,0.050000},{0.862000,0.065500,0.072500},{0.851650,0.059475,0.088875},{0.838805,0.057451,0.103744}}; alglib.mcpdcreate(3, out s); alglib.mcpdaddtrack(s, track0); alglib.mcpdaddtrack(s, track1); alglib.mcpdaddec(s, 0, 2, 0.0); alglib.mcpdaddec(s, 1, 2, 0.0); alglib.mcpdaddec(s, 2, 2, 1.0); alglib.mcpdaddec(s, 2, 0, 0.0); alglib.mcpdsolve(s); alglib.mcpdresults(s, out p, out rep); // // Hidden matrix P is equal to // ( 0.95 0.50 ) // ( 0.05 0.25 ) // ( 0.25 1.00 ) // which means that "good" loans can become past due with 5% probability, // while past due loans will become charged off with 25% probability or // return back to normal state with 50% probability. // System.Console.WriteLine("{0}", alglib.ap.format(p,2)); // EXPECTED: [[0.95,0.50,0.00],[0.05,0.25,0.00],[0.00,0.25,1.00]] System.Console.ReadLine(); return 0; }
minbleic
subpackage| minbleic_d_1 | Nonlinear optimization with bound constraints | |
| minbleic_d_2 | Nonlinear optimization with linear inequality constraints | |
| minbleic_ftrim | Nonlinear optimization by BLEIC, function with singularities | |
| minbleic_numdiff | Nonlinear optimization with bound constraints and numerical differentiation |
minbleicreport
class/************************************************************************* This structure stores optimization report: * IterationsCount number of iterations * NFEV number of gradient evaluations * TerminationType termination type (see below) TERMINATION CODES TerminationType field contains completion code, which can be: -7 gradient verification failed. See MinBLEICSetGradientCheck() for more information. -3 inconsistent constraints. Feasible point is either nonexistent or too hard to find. Try to restart optimizer with better initial approximation 1 relative function improvement is no more than EpsF. 2 relative step is no more than EpsX. 4 gradient norm is no more than EpsG 5 MaxIts steps was taken 7 stopping conditions are too stringent, further improvement is impossible, X contains best point found so far. ADDITIONAL FIELDS There are additional fields which can be used for debugging: * DebugEqErr error in the equality constraints (2-norm) * DebugFS f, calculated at projection of initial point to the feasible set * DebugFF f, calculated at the final point * DebugDX |X_start-X_final| *************************************************************************/public class minbleicreport { public int iterationscount; public int nfev; public int varidx; public int terminationtype; public double debugeqerr; public double debugfs; public double debugff; public double debugdx; public int debugfeasqpits; public int debugfeasgpaits; public int inneriterationscount; public int outeriterationscount; }
minbleicstate
class/************************************************************************* This object stores nonlinear optimizer state. You should use functions provided by MinBLEIC subpackage to work with this object *************************************************************************/public class minbleicstate { }
minbleiccreate
function/************************************************************************* BOUND CONSTRAINED OPTIMIZATION WITH ADDITIONAL LINEAR EQUALITY AND INEQUALITY CONSTRAINTS DESCRIPTION: The subroutine minimizes function F(x) of N arguments subject to any combination of: * bound constraints * linear inequality constraints * linear equality constraints REQUIREMENTS: * user must provide function value and gradient * starting point X0 must be feasible or not too far away from the feasible set * grad(f) must be Lipschitz continuous on a level set: L = { x : f(x)<=f(x0) } * function must be defined everywhere on the feasible set F USAGE: Constrained optimization if far more complex than the unconstrained one. Here we give very brief outline of the BLEIC optimizer. We strongly recommend you to read examples in the ALGLIB Reference Manual and to read ALGLIB User Guide on optimization, which is available at http://www.alglib.net/optimization/ 1. User initializes algorithm state with MinBLEICCreate() call 2. USer adds boundary and/or linear constraints by calling MinBLEICSetBC() and MinBLEICSetLC() functions. 3. User sets stopping conditions with MinBLEICSetCond(). 4. User calls MinBLEICOptimize() function which takes algorithm state and pointer (delegate, etc.) to callback function which calculates F/G. 5. User calls MinBLEICResults() to get solution 6. Optionally user may call MinBLEICRestartFrom() to solve another problem with same N but another starting point. MinBLEICRestartFrom() allows to reuse already initialized structure. INPUT PARAMETERS: N - problem dimension, N>0: * if given, only leading N elements of X are used * if not given, automatically determined from size ofX X - starting point, array[N]: * it is better to set X to a feasible point * but X can be infeasible, in which case algorithm will try to find feasible point first, using X as initial approximation. OUTPUT PARAMETERS: State - structure stores algorithm state -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleiccreate( double[] x, out minbleicstate state) public static void alglib.minbleiccreate( int n, double[] x, out minbleicstate state)
minbleiccreatef
function/************************************************************************* The subroutine is finite difference variant of MinBLEICCreate(). It uses finite differences in order to differentiate target function. Description below contains information which is specific to this function only. We recommend to read comments on MinBLEICCreate() in order to get more information about creation of BLEIC optimizer. INPUT PARAMETERS: N - problem dimension, N>0: * if given, only leading N elements of X are used * if not given, automatically determined from size of X X - starting point, array[0..N-1]. DiffStep- differentiation step, >0 OUTPUT PARAMETERS: State - structure which stores algorithm state NOTES: 1. algorithm uses 4-point central formula for differentiation. 2. differentiation step along I-th axis is equal to DiffStep*S[I] where S[] is scaling vector which can be set by MinBLEICSetScale() call. 3. we recommend you to use moderate values of differentiation step. Too large step will result in too large truncation errors, while too small step will result in too large numerical errors. 1.0E-6 can be good value to start with. 4. Numerical differentiation is very inefficient - one gradient calculation needs 4*N function evaluations. This function will work for any N - either small (1...10), moderate (10...100) or large (100...). However, performance penalty will be too severe for any N's except for small ones. We should also say that code which relies on numerical differentiation is less robust and precise. CG needs exact gradient values. Imprecise gradient may slow down convergence, especially on highly nonlinear problems. Thus we recommend to use this function for fast prototyping on small- dimensional problems only, and to implement analytical gradient as soon as possible. -- ALGLIB -- Copyright 16.05.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleiccreatef( double[] x, double diffstep, out minbleicstate state) public static void alglib.minbleiccreatef( int n, double[] x, double diffstep, out minbleicstate state)
Examples: [1]
minbleicoptimize
function/************************************************************************* This family of functions is used to launcn iterations of nonlinear optimizer These functions accept following parameters: state - algorithm state func - callback which calculates function (or merit function) value func at given point x grad - callback which calculates function (or merit function) value func and gradient grad at given point x rep - optional callback which is called after each iteration can be NULL ptr - optional pointer which is passed to func/grad/hess/jac/rep can be NULL NOTES: 1. This function has two different implementations: one which uses exact (analytical) user-supplied gradient, and one which uses function value only and numerically differentiates function in order to obtain gradient. Depending on the specific function used to create optimizer object (either MinBLEICCreate() for analytical gradient or MinBLEICCreateF() for numerical differentiation) you should choose appropriate variant of MinBLEICOptimize() - one which accepts function AND gradient or one which accepts function ONLY. Be careful to choose variant of MinBLEICOptimize() which corresponds to your optimization scheme! Table below lists different combinations of callback (function/gradient) passed to MinBLEICOptimize() and specific function used to create optimizer. | USER PASSED TO MinBLEICOptimize() CREATED WITH | function only | function and gradient ------------------------------------------------------------ MinBLEICCreateF() | work FAIL MinBLEICCreate() | FAIL work Here "FAIL" denotes inappropriate combinations of optimizer creation function and MinBLEICOptimize() version. Attemps to use such combination (for example, to create optimizer with MinBLEICCreateF() and to pass gradient information to MinCGOptimize()) will lead to exception being thrown. Either you did not pass gradient when it WAS needed or you passed gradient when it was NOT needed. -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void minbleicoptimize(minbleicstate state, ndimensional_func func, ndimensional_rep rep, object obj) public static void minbleicoptimize(minbleicstate state, ndimensional_grad grad, ndimensional_rep rep, object obj)
minbleicrestartfrom
function/************************************************************************* This subroutine restarts algorithm from new point. All optimization parameters (including constraints) are left unchanged. This function allows to solve multiple optimization problems (which must have same number of dimensions) without object reallocation penalty. INPUT PARAMETERS: State - structure previously allocated with MinBLEICCreate call. X - new starting point. -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicrestartfrom( minbleicstate state, double[] x)
minbleicresults
function/************************************************************************* BLEIC results INPUT PARAMETERS: State - algorithm state OUTPUT PARAMETERS: X - array[0..N-1], solution Rep - optimization report. You should check Rep.TerminationType in order to distinguish successful termination from unsuccessful one: * -7 gradient verification failed. See MinBLEICSetGradientCheck() for more information. * -3 inconsistent constraints. Feasible point is either nonexistent or too hard to find. Try to restart optimizer with better initial approximation * 1 relative function improvement is no more than EpsF. * 2 scaled step is no more than EpsX. * 4 scaled gradient norm is no more than EpsG. * 5 MaxIts steps was taken More information about fields of this structure can be found in the comments on MinBLEICReport datatype. -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicresults( minbleicstate state, out double[] x, out minbleicreport rep)
minbleicresultsbuf
function/************************************************************************* BLEIC results Buffered implementation of MinBLEICResults() which uses pre-allocated buffer to store X[]. If buffer size is too small, it resizes buffer. It is intended to be used in the inner cycles of performance critical algorithms where array reallocation penalty is too large to be ignored. -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicresultsbuf( minbleicstate state, ref double[] x, minbleicreport rep)
minbleicsetbc
function/************************************************************************* This function sets boundary constraints for BLEIC optimizer. Boundary constraints are inactive by default (after initial creation). They are preserved after algorithm restart with MinBLEICRestartFrom(). INPUT PARAMETERS: State - structure stores algorithm state BndL - lower bounds, array[N]. If some (all) variables are unbounded, you may specify very small number or -INF. BndU - upper bounds, array[N]. If some (all) variables are unbounded, you may specify very large number or +INF. NOTE 1: it is possible to specify BndL[i]=BndU[i]. In this case I-th variable will be "frozen" at X[i]=BndL[i]=BndU[i]. NOTE 2: this solver has following useful properties: * bound constraints are always satisfied exactly * function is evaluated only INSIDE area specified by bound constraints, even when numerical differentiation is used (algorithm adjusts nodes according to boundary constraints) -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetbc( minbleicstate state, double[] bndl, double[] bndu)
minbleicsetcond
function/************************************************************************* This function sets stopping conditions for the optimizer. INPUT PARAMETERS: State - structure which stores algorithm state EpsG - >=0 The subroutine finishes its work if the condition |v|<EpsG is satisfied, where: * |.| means Euclidian norm * v - scaled gradient vector, v[i]=g[i]*s[i] * g - gradient * s - scaling coefficients set by MinBLEICSetScale() EpsF - >=0 The subroutine finishes its work if on k+1-th iteration the condition |F(k+1)-F(k)|<=EpsF*max{|F(k)|,|F(k+1)|,1} is satisfied. EpsX - >=0 The subroutine finishes its work if on k+1-th iteration the condition |v|<=EpsX is fulfilled, where: * |.| means Euclidian norm * v - scaled step vector, v[i]=dx[i]/s[i] * dx - step vector, dx=X(k+1)-X(k) * s - scaling coefficients set by MinBLEICSetScale() MaxIts - maximum number of iterations. If MaxIts=0, the number of iterations is unlimited. Passing EpsG=0, EpsF=0 and EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic stopping criterion selection. NOTE: when SetCond() called with non-zero MaxIts, BLEIC solver may perform slightly more than MaxIts iterations. I.e., MaxIts sets non-strict limit on iterations count. -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetcond( minbleicstate state, double epsg, double epsf, double epsx, int maxits)
minbleicsetgradientcheck
function/************************************************************************* This subroutine turns on verification of the user-supplied analytic gradient: * user calls this subroutine before optimization begins * MinBLEICOptimize() is called * prior to actual optimization, for each component of parameters being optimized X[i] algorithm performs following steps: * two trial steps are made to X[i]-TestStep*S[i] and X[i]+TestStep*S[i], where X[i] is i-th component of the initial point and S[i] is a scale of i-th parameter * if needed, steps are bounded with respect to constraints on X[] * F(X) is evaluated at these trial points * we perform one more evaluation in the middle point of the interval * we build cubic model using function values and derivatives at trial points and we compare its prediction with actual value in the middle point * in case difference between prediction and actual value is higher than some predetermined threshold, algorithm stops with completion code -7; Rep.VarIdx is set to index of the parameter with incorrect derivative. * after verification is over, algorithm proceeds to the actual optimization. NOTE 1: verification needs N (parameters count) gradient evaluations. It is very costly and you should use it only for low dimensional problems, when you want to be sure that you've correctly calculated analytic derivatives. You should not use it in the production code (unless you want to check derivatives provided by some third party). NOTE 2: you should carefully choose TestStep. Value which is too large (so large that function behaviour is significantly non-cubic) will lead to false alarms. You may use different step for different parameters by means of setting scale with MinBLEICSetScale(). NOTE 3: this function may lead to false positives. In case it reports that I-th derivative was calculated incorrectly, you may decrease test step and try one more time - maybe your function changes too sharply and your step is too large for such rapidly chanding function. INPUT PARAMETERS: State - structure used to store algorithm state TestStep - verification step: * TestStep=0 turns verification off * TestStep>0 activates verification -- ALGLIB -- Copyright 15.06.2012 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetgradientcheck( minbleicstate state, double teststep)
minbleicsetlc
function/************************************************************************* This function sets linear constraints for BLEIC optimizer. Linear constraints are inactive by default (after initial creation). They are preserved after algorithm restart with MinBLEICRestartFrom(). INPUT PARAMETERS: State - structure previously allocated with MinBLEICCreate call. C - linear constraints, array[K,N+1]. Each row of C represents one constraint, either equality or inequality (see below): * first N elements correspond to coefficients, * last element corresponds to the right part. All elements of C (including right part) must be finite. CT - type of constraints, array[K]: * if CT[i]>0, then I-th constraint is C[i,*]*x >= C[i,n+1] * if CT[i]=0, then I-th constraint is C[i,*]*x = C[i,n+1] * if CT[i]<0, then I-th constraint is C[i,*]*x <= C[i,n+1] K - number of equality/inequality constraints, K>=0: * if given, only leading K elements of C/CT are used * if not given, automatically determined from sizes of C/CT NOTE 1: linear (non-bound) constraints are satisfied only approximately: * there always exists some minor violation (about Epsilon in magnitude) due to rounding errors * numerical differentiation, if used, may lead to function evaluations outside of the feasible area, because algorithm does NOT change numerical differentiation formula according to linear constraints. If you want constraints to be satisfied exactly, try to reformulate your problem in such manner that all constraints will become boundary ones (this kind of constraints is always satisfied exactly, both in the final solution and in all intermediate points). -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetlc( minbleicstate state, double[,] c, int[] ct) public static void alglib.minbleicsetlc( minbleicstate state, double[,] c, int[] ct, int k)
Examples: [1]
minbleicsetprecdefault
function/************************************************************************* Modification of the preconditioner: preconditioning is turned off. INPUT PARAMETERS: State - structure which stores algorithm state -- ALGLIB -- Copyright 13.10.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetprecdefault(minbleicstate state)
minbleicsetprecdiag
function/************************************************************************* Modification of the preconditioner: diagonal of approximate Hessian is used. INPUT PARAMETERS: State - structure which stores algorithm state D - diagonal of the approximate Hessian, array[0..N-1], (if larger, only leading N elements are used). NOTE 1: D[i] should be positive. Exception will be thrown otherwise. NOTE 2: you should pass diagonal of approximate Hessian - NOT ITS INVERSE. -- ALGLIB -- Copyright 13.10.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetprecdiag( minbleicstate state, double[] d)
minbleicsetprecscale
function/************************************************************************* Modification of the preconditioner: scale-based diagonal preconditioning. This preconditioning mode can be useful when you don't have approximate diagonal of Hessian, but you know that your variables are badly scaled (for example, one variable is in [1,10], and another in [1000,100000]), and most part of the ill-conditioning comes from different scales of vars. In this case simple scale-based preconditioner, with H[i] = 1/(s[i]^2), can greatly improve convergence. IMPRTANT: you should set scale of your variables with MinBLEICSetScale() call (before or after MinBLEICSetPrecScale() call). Without knowledge of the scale of your variables scale-based preconditioner will be just unit matrix. INPUT PARAMETERS: State - structure which stores algorithm state -- ALGLIB -- Copyright 13.10.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetprecscale(minbleicstate state)
minbleicsetscale
function/************************************************************************* This function sets scaling coefficients for BLEIC optimizer. ALGLIB optimizers use scaling matrices to test stopping conditions (step size and gradient are scaled before comparison with tolerances). Scale of the I-th variable is a translation invariant measure of: a) "how large" the variable is b) how large the step should be to make significant changes in the function Scaling is also used by finite difference variant of the optimizer - step along I-th axis is equal to DiffStep*S[I]. In most optimizers (and in the BLEIC too) scaling is NOT a form of preconditioning. It just affects stopping conditions. You should set preconditioner by separate call to one of the MinBLEICSetPrec...() functions. There is a special preconditioning mode, however, which uses scaling coefficients to form diagonal preconditioning matrix. You can turn this mode on, if you want. But you should understand that scaling is not the same thing as preconditioning - these are two different, although related forms of tuning solver. INPUT PARAMETERS: State - structure stores algorithm state S - array[N], non-zero scaling coefficients S[i] may be negative, sign doesn't matter. -- ALGLIB -- Copyright 14.01.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetscale( minbleicstate state, double[] s)
minbleicsetstpmax
function/************************************************************************* This function sets maximum step length IMPORTANT: this feature is hard to combine with preconditioning. You can't set upper limit on step length, when you solve optimization problem with linear (non-boundary) constraints AND preconditioner turned on. When non-boundary constraints are present, you have to either a) use preconditioner, or b) use upper limit on step length. YOU CAN'T USE BOTH! In this case algorithm will terminate with appropriate error code. INPUT PARAMETERS: State - structure which stores algorithm state StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't want to limit step length. Use this subroutine when you optimize target function which contains exp() or other fast growing functions, and optimization algorithm makes too large steps which lead to overflow. This function allows us to reject steps that are too large (and therefore expose us to the possible overflow) without actually calculating function value at the x+stp*d. -- ALGLIB -- Copyright 02.04.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetstpmax( minbleicstate state, double stpmax)
minbleicsetxrep
function/************************************************************************* This function turns on/off reporting. INPUT PARAMETERS: State - structure which stores algorithm state NeedXRep- whether iteration reports are needed or not If NeedXRep is True, algorithm will call rep() callback function if it is provided to MinBLEICOptimize(). -- ALGLIB -- Copyright 28.11.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.minbleicsetxrep( minbleicstate state, bool needxrep)
public static void function1_grad(double[] x, ref double func, double[] grad, object obj) { // this callback calculates f(x0,x1) = 100*(x0+3)^4 + (x1-3)^4 // and its derivatives df/d0 and df/dx1 func = 100*System.Math.Pow(x[0]+3,4) + System.Math.Pow(x[1]-3,4); grad[0] = 400*System.Math.Pow(x[0]+3,3); grad[1] = 4*System.Math.Pow(x[1]-3,3); } public static int Main(string[] args) { // // This example demonstrates minimization of f(x,y) = 100*(x+3)^4+(y-3)^4 // subject to bound constraints -1<=x<=+1, -1<=y<=+1, using BLEIC optimizer. // double[] x = new double[]{0,0}; double[] bndl = new double[]{-1,-1}; double[] bndu = new double[]{+1,+1}; alglib.minbleicstate state; alglib.minbleicreport rep; // // These variables define stopping conditions for the optimizer. // // We use very simple condition - |g|<=epsg // double epsg = 0.000001; double epsf = 0; double epsx = 0; int maxits = 0; // // Now we are ready to actually optimize something: // * first we create optimizer // * we add boundary constraints // * we tune stopping conditions // * and, finally, optimize and obtain results... // alglib.minbleiccreate(x, out state); alglib.minbleicsetbc(state, bndl, bndu); alglib.minbleicsetcond(state, epsg, epsf, epsx, maxits); alglib.minbleicoptimize(state, function1_grad, null, null); alglib.minbleicresults(state, out x, out rep); // // ...and evaluate these results // System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4 System.Console.WriteLine("{0}", alglib.ap.format(x,2)); // EXPECTED: [-1,1] System.Console.ReadLine(); return 0; }
public static void function1_grad(double[] x, ref double func, double[] grad, object obj) { // this callback calculates f(x0,x1) = 100*(x0+3)^4 + (x1-3)^4 // and its derivatives df/d0 and df/dx1 func = 100*System.Math.Pow(x[0]+3,4) + System.Math.Pow(x[1]-3,4); grad[0] = 400*System.Math.Pow(x[0]+3,3); grad[1] = 4*System.Math.Pow(x[1]-3,3); } public static int Main(string[] args) { // // This example demonstrates minimization of f(x,y) = 100*(x+3)^4+(y-3)^4 // subject to inequality constraints: // * x>=2 (posed as general linear constraint), // * x+y>=6 // using BLEIC optimizer. // double[] x = new double[]{5,5}; double[,] c = new double[,]{{1,0,2},{1,1,6}}; int[] ct = new int[]{1,1}; alglib.minbleicstate state; alglib.minbleicreport rep; // // These variables define stopping conditions for the optimizer. // // We use very simple condition - |g|<=epsg // double epsg = 0.000001; double epsf = 0; double epsx = 0; int maxits = 0; // // Now we are ready to actually optimize something: // * first we create optimizer // * we add linear constraints // * we tune stopping conditions // * and, finally, optimize and obtain results... // alglib.minbleiccreate(x, out state); alglib.minbleicsetlc(state, c, ct); alglib.minbleicsetcond(state, epsg, epsf, epsx, maxits); alglib.minbleicoptimize(state, function1_grad, null, null); alglib.minbleicresults(state, out x, out rep); // // ...and evaluate these results // System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4 System.Console.WriteLine("{0}", alglib.ap.format(x,2)); // EXPECTED: [2,4] System.Console.ReadLine(); return 0; }
public static void s1_grad(double[] x, ref double func, double[] grad, object obj) { // // this callback calculates f(x) = (1+x)^(-0.2) + (1-x)^(-0.3) + 1000*x and its gradient. // // function is trimmed when we calculate it near the singular points or outside of the [-1,+1]. // Note that we do NOT calculate gradient in this case. // if( (x[0]<=-0.999999999999) || (x[0]>=+0.999999999999) ) { func = 1.0E+300; return; } func = System.Math.Pow(1+x[0],-0.2) + System.Math.Pow(1-x[0],-0.3) + 1000*x[0]; grad[0] = -0.2*System.Math.Pow(1+x[0],-1.2) +0.3*System.Math.Pow(1-x[0],-1.3) + 1000; } public static int Main(string[] args) { // // This example demonstrates minimization of f(x) = (1+x)^(-0.2) + (1-x)^(-0.3) + 1000*x. // // This function is undefined outside of (-1,+1) and has singularities at x=-1 and x=+1. // Special technique called "function trimming" allows us to solve this optimization problem // - without using boundary constraints! // // See http://www.alglib.net/optimization/tipsandtricks.php#ftrimming for more information // on this subject. // double[] x = new double[]{0}; double epsg = 1.0e-6; double epsf = 0; double epsx = 0; int maxits = 0; alglib.minbleicstate state; alglib.minbleicreport rep; alglib.minbleiccreate(x, out state); alglib.minbleicsetcond(state, epsg, epsf, epsx, maxits); alglib.minbleicoptimize(state, s1_grad, null, null); alglib.minbleicresults(state, out x, out rep); System.Console.WriteLine("{0}", alglib.ap.format(x,5)); // EXPECTED: [-0.99917305] System.Console.ReadLine(); return 0; }
public static void function1_func(double[] x, ref double func, object obj) { // this callback calculates f(x0,x1) = 100*(x0+3)^4 + (x1-3)^4 func = 100*System.Math.Pow(x[0]+3,4) + System.Math.Pow(x[1]-3,4); } public static int Main(string[] args) { // // This example demonstrates minimization of f(x,y) = 100*(x+3)^4+(y-3)^4 // subject to bound constraints -1<=x<=+1, -1<=y<=+1, using BLEIC optimizer. // double[] x = new double[]{0,0}; double[] bndl = new double[]{-1,-1}; double[] bndu = new double[]{+1,+1}; alglib.minbleicstate state; alglib.minbleicreport rep; // // These variables define stopping conditions for the optimizer. // // We use very simple condition - |g|<=epsg // double epsg = 0.000001; double epsf = 0; double epsx = 0; int maxits = 0; // // This variable contains differentiation step // double diffstep = 1.0e-6; // // Now we are ready to actually optimize something: // * first we create optimizer // * we add boundary constraints // * we tune stopping conditions // * and, finally, optimize and obtain results... // alglib.minbleiccreatef(x, diffstep, out state); alglib.minbleicsetbc(state, bndl, bndu); alglib.minbleicsetcond(state, epsg, epsf, epsx, maxits); alglib.minbleicoptimize(state, function1_func, null, null); alglib.minbleicresults(state, out x, out rep); // // ...and evaluate these results // System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4 System.Console.WriteLine("{0}", alglib.ap.format(x,2)); // EXPECTED: [-1,1] System.Console.ReadLine(); return 0; }
mincg
subpackage| mincg_d_1 | Nonlinear optimization by CG | |
| mincg_d_2 | Nonlinear optimization with additional settings and restarts | |
| mincg_ftrim | Nonlinear optimization by CG, function with singularities | |
| mincg_numdiff | Nonlinear optimization by CG with numerical differentiation |
mincgreport
class/************************************************************************* *************************************************************************/public class mincgreport { public int iterationscount; public int nfev; public int varidx; public int terminationtype; }
mincgstate
class/************************************************************************* This object stores state of the nonlinear CG optimizer. You should use ALGLIB functions to work with this object. *************************************************************************/public class mincgstate { }
mincgcreate
function/************************************************************************* NONLINEAR CONJUGATE GRADIENT METHOD DESCRIPTION: The subroutine minimizes function F(x) of N arguments by using one of the nonlinear conjugate gradient methods. These CG methods are globally convergent (even on non-convex functions) as long as grad(f) is Lipschitz continuous in a some neighborhood of the L = { x : f(x)<=f(x0) }. REQUIREMENTS: Algorithm will request following information during its operation: * function value F and its gradient G (simultaneously) at given point X USAGE: 1. User initializes algorithm state with MinCGCreate() call 2. User tunes solver parameters with MinCGSetCond(), MinCGSetStpMax() and other functions 3. User calls MinCGOptimize() function which takes algorithm state and pointer (delegate, etc.) to callback function which calculates F/G. 4. User calls MinCGResults() to get solution 5. Optionally, user may call MinCGRestartFrom() to solve another problem with same N but another starting point and/or another function. MinCGRestartFrom() allows to reuse already initialized structure. INPUT PARAMETERS: N - problem dimension, N>0: * if given, only leading N elements of X are used * if not given, automatically determined from size of X X - starting point, array[0..N-1]. OUTPUT PARAMETERS: State - structure which stores algorithm state -- ALGLIB -- Copyright 25.03.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgcreate(double[] x, out mincgstate state) public static void alglib.mincgcreate( int n, double[] x, out mincgstate state)
mincgcreatef
function/************************************************************************* The subroutine is finite difference variant of MinCGCreate(). It uses finite differences in order to differentiate target function. Description below contains information which is specific to this function only. We recommend to read comments on MinCGCreate() in order to get more information about creation of CG optimizer. INPUT PARAMETERS: N - problem dimension, N>0: * if given, only leading N elements of X are used * if not given, automatically determined from size of X X - starting point, array[0..N-1]. DiffStep- differentiation step, >0 OUTPUT PARAMETERS: State - structure which stores algorithm state NOTES: 1. algorithm uses 4-point central formula for differentiation. 2. differentiation step along I-th axis is equal to DiffStep*S[I] where S[] is scaling vector which can be set by MinCGSetScale() call. 3. we recommend you to use moderate values of differentiation step. Too large step will result in too large truncation errors, while too small step will result in too large numerical errors. 1.0E-6 can be good value to start with. 4. Numerical differentiation is very inefficient - one gradient calculation needs 4*N function evaluations. This function will work for any N - either small (1...10), moderate (10...100) or large (100...). However, performance penalty will be too severe for any N's except for small ones. We should also say that code which relies on numerical differentiation is less robust and precise. L-BFGS needs exact gradient values. Imprecise gradient may slow down convergence, especially on highly nonlinear problems. Thus we recommend to use this function for fast prototyping on small- dimensional problems only, and to implement analytical gradient as soon as possible. -- ALGLIB -- Copyright 16.05.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgcreatef( double[] x, double diffstep, out mincgstate state) public static void alglib.mincgcreatef( int n, double[] x, double diffstep, out mincgstate state)
Examples: [1]
mincgoptimize
function/************************************************************************* This family of functions is used to launcn iterations of nonlinear optimizer These functions accept following parameters: state - algorithm state func - callback which calculates function (or merit function) value func at given point x grad - callback which calculates function (or merit function) value func and gradient grad at given point x rep - optional callback which is called after each iteration can be NULL ptr - optional pointer which is passed to func/grad/hess/jac/rep can be NULL NOTES: 1. This function has two different implementations: one which uses exact (analytical) user-supplied gradient, and one which uses function value only and numerically differentiates function in order to obtain gradient. Depending on the specific function used to create optimizer object (either MinCGCreate() for analytical gradient or MinCGCreateF() for numerical differentiation) you should choose appropriate variant of MinCGOptimize() - one which accepts function AND gradient or one which accepts function ONLY. Be careful to choose variant of MinCGOptimize() which corresponds to your optimization scheme! Table below lists different combinations of callback (function/gradient) passed to MinCGOptimize() and specific function used to create optimizer. | USER PASSED TO MinCGOptimize() CREATED WITH | function only | function and gradient ------------------------------------------------------------ MinCGCreateF() | work FAIL MinCGCreate() | FAIL work Here "FAIL" denotes inappropriate combinations of optimizer creation function and MinCGOptimize() version. Attemps to use such combination (for example, to create optimizer with MinCGCreateF() and to pass gradient information to MinCGOptimize()) will lead to exception being thrown. Either you did not pass gradient when it WAS needed or you passed gradient when it was NOT needed. -- ALGLIB -- Copyright 20.04.2009 by Bochkanov Sergey *************************************************************************/public static void mincgoptimize(mincgstate state, ndimensional_func func, ndimensional_rep rep, object obj) public static void mincgoptimize(mincgstate state, ndimensional_grad grad, ndimensional_rep rep, object obj)
mincgrestartfrom
function/************************************************************************* This subroutine restarts CG algorithm from new point. All optimization parameters are left unchanged. This function allows to solve multiple optimization problems (which must have same number of dimensions) without object reallocation penalty. INPUT PARAMETERS: State - structure used to store algorithm state. X - new starting point. -- ALGLIB -- Copyright 30.07.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgrestartfrom(mincgstate state, double[] x)
Examples: [1]
mincgresults
function/************************************************************************* Conjugate gradient results INPUT PARAMETERS: State - algorithm state OUTPUT PARAMETERS: X - array[0..N-1], solution Rep - optimization report: * Rep.TerminationType completetion code: * -7 gradient verification failed. See MinCGSetGradientCheck() for more information. * 1 relative function improvement is no more than EpsF. * 2 relative step is no more than EpsX. * 4 gradient norm is no more than EpsG * 5 MaxIts steps was taken * 7 stopping conditions are too stringent, further improvement is impossible, we return best X found so far * 8 terminated by user * Rep.IterationsCount contains iterations count * NFEV countains number of function calculations -- ALGLIB -- Copyright 20.04.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgresults( mincgstate state, out double[] x, out mincgreport rep)
mincgresultsbuf
function/************************************************************************* Conjugate gradient results Buffered implementation of MinCGResults(), which uses pre-allocated buffer to store X[]. If buffer size is too small, it resizes buffer. It is intended to be used in the inner cycles of performance critical algorithms where array reallocation penalty is too large to be ignored. -- ALGLIB -- Copyright 20.04.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgresultsbuf( mincgstate state, ref double[] x, mincgreport rep)
mincgsetcgtype
function/************************************************************************* This function sets CG algorithm. INPUT PARAMETERS: State - structure which stores algorithm state CGType - algorithm type: * -1 automatic selection of the best algorithm * 0 DY (Dai and Yuan) algorithm * 1 Hybrid DY-HS algorithm -- ALGLIB -- Copyright 02.04.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgsetcgtype(mincgstate state, int cgtype)
mincgsetcond
function/************************************************************************* This function sets stopping conditions for CG optimization algorithm. INPUT PARAMETERS: State - structure which stores algorithm state EpsG - >=0 The subroutine finishes its work if the condition |v|<EpsG is satisfied, where: * |.| means Euclidian norm * v - scaled gradient vector, v[i]=g[i]*s[i] * g - gradient * s - scaling coefficients set by MinCGSetScale() EpsF - >=0 The subroutine finishes its work if on k+1-th iteration the condition |F(k+1)-F(k)|<=EpsF*max{|F(k)|,|F(k+1)|,1} is satisfied. EpsX - >=0 The subroutine finishes its work if on k+1-th iteration the condition |v|<=EpsX is fulfilled, where: * |.| means Euclidian norm * v - scaled step vector, v[i]=dx[i]/s[i] * dx - ste pvector, dx=X(k+1)-X(k) * s - scaling coefficients set by MinCGSetScale() MaxIts - maximum number of iterations. If MaxIts=0, the number of iterations is unlimited. Passing EpsG=0, EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic stopping criterion selection (small EpsX). -- ALGLIB -- Copyright 02.04.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgsetcond( mincgstate state, double epsg, double epsf, double epsx, int maxits)
mincgsetgradientcheck
function/************************************************************************* This subroutine turns on verification of the user-supplied analytic gradient: * user calls this subroutine before optimization begins * MinCGOptimize() is called * prior to actual optimization, for each component of parameters being optimized X[i] algorithm performs following steps: * two trial steps are made to X[i]-TestStep*S[i] and X[i]+TestStep*S[i], where X[i] is i-th component of the initial point and S[i] is a scale of i-th parameter * F(X) is evaluated at these trial points * we perform one more evaluation in the middle point of the interval * we build cubic model using function values and derivatives at trial points and we compare its prediction with actual value in the middle point * in case difference between prediction and actual value is higher than some predetermined threshold, algorithm stops with completion code -7; Rep.VarIdx is set to index of the parameter with incorrect derivative. * after verification is over, algorithm proceeds to the actual optimization. NOTE 1: verification needs N (parameters count) gradient evaluations. It is very costly and you should use it only for low dimensional problems, when you want to be sure that you've correctly calculated analytic derivatives. You should not use it in the production code (unless you want to check derivatives provided by some third party). NOTE 2: you should carefully choose TestStep. Value which is too large (so large that function behaviour is significantly non-cubic) will lead to false alarms. You may use different step for different parameters by means of setting scale with MinCGSetScale(). NOTE 3: this function may lead to false positives. In case it reports that I-th derivative was calculated incorrectly, you may decrease test step and try one more time - maybe your function changes too sharply and your step is too large for such rapidly chanding function. INPUT PARAMETERS: State - structure used to store algorithm state TestStep - verification step: * TestStep=0 turns verification off * TestStep>0 activates verification -- ALGLIB -- Copyright 31.05.2012 by Bochkanov Sergey *************************************************************************/public static void alglib.mincgsetgradientcheck( mincgstate state, double teststep)
mincgsetprecdefault