This is octave.info, produced by makeinfo version 5.2 from octave.texi. START-INFO-DIR-ENTRY * Octave: (octave). Interactive language for numerical computations. END-INFO-DIR-ENTRY Copyright (C) 1996-2013 John W. Eaton. Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions.  File: octave.info, Node: Explicit and Implicit Conversions, Prev: Creating Permutation Matrices, Up: Basic Usage 21.1.3 Explicit and Implicit Conversions ---------------------------------------- The diagonal and permutation matrices are special objects in their own right. A number of operations and built-in functions are defined for these matrices to use special, more efficient code than would be used for a full matrix in the same place. Examples are given in further sections. To facilitate smooth mixing with full matrices, backward compatibility, and compatibility with MATLAB, the diagonal and permutation matrices should allow any operation that works on full matrices, and will either treat it specially, or implicitly convert themselves to full matrices. Instances include matrix indexing, except for extracting a single element or a leading submatrix, indexed assignment, or applying most mapper functions, such as "exp". An explicit conversion to a full matrix can be requested using the built-in function "full". It should also be noted that the diagonal and permutation matrix objects will cache the result of the conversion after it is first requested (explicitly or implicitly), so that subsequent conversions will be very cheap.  File: octave.info, Node: Matrix Algebra, Next: Function Support, Prev: Basic Usage, Up: Diagonal and Permutation Matrices 21.2 Linear Algebra with Diagonal/Permutation Matrices ====================================================== As has been already said, diagonal and permutation matrices make it possible to use efficient algorithms while preserving natural linear algebra syntax. This section describes in detail the operations that are treated specially when performed on these special matrix objects. * Menu: * Expressions Involving Diagonal Matrices:: * Expressions Involving Permutation Matrices::  File: octave.info, Node: Expressions Involving Diagonal Matrices, Next: Expressions Involving Permutation Matrices, Up: Matrix Algebra 21.2.1 Expressions Involving Diagonal Matrices ---------------------------------------------- Assume D is a diagonal matrix. If M is a full matrix, then 'D*M' will scale the rows of M. That means, if 'S = D*M', then for each pair of indices i,j it holds S(i,j) = D(i,i) * M(i,j). Similarly, 'M*D' will do a column scaling. The matrix D may also be rectangular, m-by-n where 'm != n'. If 'm < n', then the expression 'D*M' is equivalent to D(:,1:m) * M(1:m,:), i.e., trailing 'n-m' rows of M are ignored. If 'm > n', then 'D*M' is equivalent to [D(1:n,n) * M; zeros(m-n, columns (M))], i.e., null rows are appended to the result. The situation for right-multiplication 'M*D' is analogous. The expressions 'D \ M' and 'M / D' perform inverse scaling. They are equivalent to solving a diagonal (or rectangular diagonal) in a least-squares minimum-norm sense. In exact arithmetic, this is equivalent to multiplying by a pseudoinverse. The pseudoinverse of a rectangular diagonal matrix is again a rectangular diagonal matrix with swapped dimensions, where each nonzero diagonal element is replaced by its reciprocal. The matrix division algorithms do, in fact, use division rather than multiplication by reciprocals for better numerical accuracy; otherwise, they honor the above definition. Note that a diagonal matrix is never truncated due to ill-conditioning; otherwise, it would not be of much use for scaling. This is typically consistent with linear algebra needs. A full matrix that only happens to be diagonal (and is thus not a special object) is of course treated normally. Multiplication and division by diagonal matrices work efficiently also when combined with sparse matrices, i.e., 'D*S', where D is a diagonal matrix and S is a sparse matrix scales the rows of the sparse matrix and returns a sparse matrix. The expressions 'S*D', 'D\S', 'S/D' work analogically. If D1 and D2 are both diagonal matrices, then the expressions D1 + D2 D1 - D2 D1 * D2 D1 / D2 D1 \ D2 again produce diagonal matrices, provided that normal dimension matching rules are obeyed. The relations used are same as described above. Also, a diagonal matrix D can be multiplied or divided by a scalar, or raised to a scalar power if it is square, producing diagonal matrix result in all cases. A diagonal matrix can also be transposed or conjugate-transposed, giving the expected result. Extracting a leading submatrix of a diagonal matrix, i.e., 'D(1:m,1:n)', will produce a diagonal matrix, other indexing expressions will implicitly convert to full matrix. Adding a diagonal matrix to a full matrix only operates on the diagonal elements. Thus, A = A + eps * eye (n) is an efficient method of augmenting the diagonal of a matrix. Subtraction works analogically. When involved in expressions with other element-by-element operators, '.*', './', '.\' or '.^', an implicit conversion to full matrix will take place. This is not always strictly necessary but chosen to facilitate better consistency with MATLAB.  File: octave.info, Node: Expressions Involving Permutation Matrices, Prev: Expressions Involving Diagonal Matrices, Up: Matrix Algebra 21.2.2 Expressions Involving Permutation Matrices ------------------------------------------------- If P is a permutation matrix and M a matrix, the expression 'P*M' will permute the rows of M. Similarly, 'M*P' will yield a column permutation. Matrix division 'P\M' and 'M/P' can be used to do inverse permutation. The previously described syntax for creating permutation matrices can actually help an user to understand the connection between a permutation matrix and a permuting vector. Namely, the following holds, where 'I = eye (n)' is an identity matrix: I(p,:) * M = (I*M) (p,:) = M(p,:) Similarly, M * I(:,p) = (M*I) (:,p) = M(:,p) The expressions 'I(p,:)' and 'I(:,p)' are permutation matrices. A permutation matrix can be transposed (or conjugate-transposed, which is the same, because a permutation matrix is never complex), inverting the permutation, or equivalently, turning a row-permutation matrix into a column-permutation one. For permutation matrices, transpose is equivalent to inversion, thus 'P\M' is equivalent to 'P'*M'. Transpose of a permutation matrix (or inverse) is a constant-time operation, flipping only a flag internally, and thus the choice between the two above equivalent expressions for inverse permuting is completely up to the user's taste. Multiplication and division by permutation matrices works efficiently also when combined with sparse matrices, i.e., 'P*S', where P is a permutation matrix and S is a sparse matrix permutes the rows of the sparse matrix and returns a sparse matrix. The expressions 'S*P', 'P\S', 'S/P' work analogically. Two permutation matrices can be multiplied or divided (if their sizes match), performing a composition of permutations. Also a permutation matrix can be indexed by a permutation vector (or two vectors), giving again a permutation matrix. Any other operations do not generally yield a permutation matrix and will thus trigger the implicit conversion.  File: octave.info, Node: Function Support, Next: Example Code, Prev: Matrix Algebra, Up: Diagonal and Permutation Matrices 21.3 Functions That Are Aware of These Matrices =============================================== This section lists the built-in functions that are aware of diagonal and permutation matrices on input, or can return them as output. Passed to other functions, these matrices will in general trigger an implicit conversion. (Of course, user-defined dynamically linked functions may also work with diagonal or permutation matrices). * Menu: * Diagonal Matrix Functions:: * Permutation Matrix Functions::  File: octave.info, Node: Diagonal Matrix Functions, Next: Permutation Matrix Functions, Up: Function Support 21.3.1 Diagonal Matrix Functions -------------------------------- "inv" and "pinv" can be applied to a diagonal matrix, yielding again a diagonal matrix. "det" will use an efficient straightforward calculation when given a diagonal matrix, as well as "cond". The following mapper functions can be applied to a diagonal matrix without converting it to a full one: "abs", "real", "imag", "conj", "sqrt". A diagonal matrix can also be returned from the "balance" and "svd" functions. The "sparse" function will convert a diagonal matrix efficiently to a sparse matrix.  File: octave.info, Node: Permutation Matrix Functions, Prev: Diagonal Matrix Functions, Up: Function Support 21.3.2 Permutation Matrix Functions ----------------------------------- "inv" and "pinv" will invert a permutation matrix, preserving its specialness. "det" can be applied to a permutation matrix, efficiently calculating the sign of the permutation (which is equal to the determinant). A permutation matrix can also be returned from the built-in functions "lu" and "qr", if a pivoted factorization is requested. The "sparse" function will convert a permutation matrix efficiently to a sparse matrix. The "find" function will also work efficiently with a permutation matrix, making it possible to conveniently obtain the permutation indices.  File: octave.info, Node: Example Code, Next: Zeros Treatment, Prev: Function Support, Up: Diagonal and Permutation Matrices 21.4 Examples of Usage ====================== The following can be used to solve a linear system 'A*x = b' using the pivoted LU factorization: [L, U, P] = lu (A); ## now L*U = P*A x = U \ L \ P*b; This is one way to normalize columns of a matrix X to unit norm: s = norm (X, "columns"); X /= diag (s); The same can also be accomplished with broadcasting (*note Broadcasting::): s = norm (X, "columns"); X ./= s; The following expression is a way to efficiently calculate the sign of a permutation, given by a permutation vector P. It will also work in earlier versions of Octave, but slowly. det (eye (length (p))(p, :)) Finally, here's how to solve a linear system 'A*x = b' with Tikhonov regularization (ridge regression) using SVD (a skeleton only): m = rows (A); n = columns (A); [U, S, V] = svd (A); ## determine the regularization factor alpha ## alpha = ... ## transform to orthogonal basis b = U'*b; ## Use the standard formula, replacing A with S. ## S is diagonal, so the following will be very fast and accurate. x = (S'*S + alpha^2 * eye (n)) \ (S' * b); ## transform to solution basis x = V*x;  File: octave.info, Node: Zeros Treatment, Prev: Example Code, Up: Diagonal and Permutation Matrices 21.5 Differences in Treatment of Zero Elements ============================================== Making diagonal and permutation matrices special matrix objects in their own right and the consequent usage of smarter algorithms for certain operations implies, as a side effect, small differences in treating zeros. The contents of this section apply also to sparse matrices, discussed in the following chapter. (*note Sparse Matrices::) The IEEE floating point standard defines the result of the expressions '0*Inf' and '0*NaN' as 'NaN'. This is widely agreed to be a good compromise. Numerical software dealing with structured and sparse matrices (including Octave) however, almost always makes a distinction between a "numerical zero" and an "assumed zero". A "numerical zero" is a zero value occurring in a place where any floating-point value could occur. It is normally stored somewhere in memory as an explicit value. An "assumed zero", on the contrary, is a zero matrix element implied by the matrix structure (diagonal, triangular) or a sparsity pattern; its value is usually not stored explicitly anywhere, but is implied by the underlying data structure. The primary distinction is that an assumed zero, when multiplied by any number, or divided by any nonzero number, yields *always* a zero, even when, e.g., multiplied by 'Inf' or divided by 'NaN'. The reason for this behavior is that the numerical multiplication is not actually performed anywhere by the underlying algorithm; the result is just assumed to be zero. Equivalently, one can say that the part of the computation involving assumed zeros is performed symbolically, not numerically. This behavior not only facilitates the most straightforward and efficient implementation of algorithms, but also preserves certain useful invariants, like: * scalar * diagonal matrix is a diagonal matrix * sparse matrix / scalar preserves the sparsity pattern * permutation matrix * matrix is equivalent to permuting rows all of these natural mathematical truths would be invalidated by treating assumed zeros as numerical ones. Note that MATLAB does not strictly follow this principle and converts assumed zeros to numerical zeros in certain cases, while not doing so in other cases. As of today, there are no intentions to mimic such behavior in Octave. Examples of effects of assumed zeros vs. numerical zeros: Inf * eye (3) => Inf 0 0 0 Inf 0 0 0 Inf Inf * speye (3) => Compressed Column Sparse (rows = 3, cols = 3, nnz = 3 [33%]) (1, 1) -> Inf (2, 2) -> Inf (3, 3) -> Inf Inf * full (eye (3)) => Inf NaN NaN NaN Inf NaN NaN NaN Inf diag (1:3) * [NaN; 1; 1] => NaN 2 3 sparse (1:3,1:3,1:3) * [NaN; 1; 1] => NaN 2 3 [1,0,0;0,2,0;0,0,3] * [NaN; 1; 1] => NaN NaN NaN  File: octave.info, Node: Sparse Matrices, Next: Numerical Integration, Prev: Diagonal and Permutation Matrices, Up: Top 22 Sparse Matrices ****************** * Menu: * Basics:: Creation and Manipulation of Sparse Matrices * Sparse Linear Algebra:: Linear Algebra on Sparse Matrices * Iterative Techniques:: Iterative Techniques * Real Life Example:: Using Sparse Matrices  File: octave.info, Node: Basics, Next: Sparse Linear Algebra, Up: Sparse Matrices 22.1 Creation and Manipulation of Sparse Matrices ================================================= The size of mathematical problems that can be treated at any particular time is generally limited by the available computing resources. Both, the speed of the computer and its available memory place limitation on the problem size. There are many classes of mathematical problems which give rise to matrices, where a large number of the elements are zero. In this case it makes sense to have a special matrix type to handle this class of problems where only the nonzero elements of the matrix are stored. Not only does this reduce the amount of memory to store the matrix, but it also means that operations on this type of matrix can take advantage of the a priori knowledge of the positions of the nonzero elements to accelerate their calculations. A matrix type that stores only the nonzero elements is generally called sparse. It is the purpose of this document to discuss the basics of the storage and creation of sparse matrices and the fundamental operations on them. * Menu: * Storage of Sparse Matrices:: * Creating Sparse Matrices:: * Information:: * Operators and Functions::  File: octave.info, Node: Storage of Sparse Matrices, Next: Creating Sparse Matrices, Up: Basics 22.1.1 Storage of Sparse Matrices --------------------------------- It is not strictly speaking necessary for the user to understand how sparse matrices are stored. However, such an understanding will help to get an understanding of the size of sparse matrices. Understanding the storage technique is also necessary for those users wishing to create their own oct-files. There are many different means of storing sparse matrix data. What all of the methods have in common is that they attempt to reduce the complexity and storage given a priori knowledge of the particular class of problems that will be solved. A good summary of the available techniques for storing sparse matrix is given by Saad (1). With full matrices, knowledge of the point of an element of the matrix within the matrix is implied by its position in the computers memory. However, this is not the case for sparse matrices, and so the positions of the nonzero elements of the matrix must equally be stored. An obvious way to do this is by storing the elements of the matrix as triplets, with two elements being their position in the array (rows and column) and the third being the data itself. This is conceptually easy to grasp, but requires more storage than is strictly needed. The storage technique used within Octave is the compressed column format. It is similar to the Yale format. (2) In this format the position of each element in a row and the data are stored as previously. However, if we assume that all elements in the same column are stored adjacent in the computers memory, then we only need to store information on the number of nonzero elements in each column, rather than their positions. Thus assuming that the matrix has more nonzero elements than there are columns in the matrix, we win in terms of the amount of memory used. In fact, the column index contains one more element than the number of columns, with the first element always being zero. The advantage of this is a simplification in the code, in that there is no special case for the first or last columns. A short example, demonstrating this in C is. for (j = 0; j < nc; j++) for (i = cidx(j); i < cidx(j+1); i++) printf ("nonzero element (%i,%i) is %d\n", ridx(i), j, data(i)); A clear understanding might be had by considering an example of how the above applies to an example matrix. Consider the matrix 1 2 0 0 0 0 0 3 0 0 0 4 The nonzero elements of this matrix are (1, 1) => 1 (1, 2) => 2 (2, 4) => 3 (3, 4) => 4 This will be stored as three vectors CIDX, RIDX and DATA, representing the column indexing, row indexing and data respectively. The contents of these three vectors for the above matrix will be CIDX = [0, 1, 2, 2, 4] RIDX = [0, 0, 1, 2] DATA = [1, 2, 3, 4] Note that this is the representation of these elements with the first row and column assumed to start at zero, while in Octave itself the row and column indexing starts at one. Thus the number of elements in the I-th column is given by 'CIDX (I + 1) - CIDX (I)'. Although Octave uses a compressed column format, it should be noted that compressed row formats are equally possible. However, in the context of mixed operations between mixed sparse and dense matrices, it makes sense that the elements of the sparse matrices are in the same order as the dense matrices. Octave stores dense matrices in column major ordering, and so sparse matrices are equally stored in this manner. A further constraint on the sparse matrix storage used by Octave is that all elements in the rows are stored in increasing order of their row index, which makes certain operations faster. However, it imposes the need to sort the elements on the creation of sparse matrices. Having disordered elements is potentially an advantage in that it makes operations such as concatenating two sparse matrices together easier and faster, however it adds complexity and speed problems elsewhere. ---------- Footnotes ---------- (1) Y. Saad "SPARSKIT: A basic toolkit for sparse matrix computation", 1994, (2)  File: octave.info, Node: Creating Sparse Matrices, Next: Information, Prev: Storage of Sparse Matrices, Up: Basics 22.1.2 Creating Sparse Matrices ------------------------------- There are several means to create sparse matrix. Returned from a function There are many functions that directly return sparse matrices. These include "speye", "sprand", "diag", etc. Constructed from matrices or vectors The function "sparse" allows a sparse matrix to be constructed from three vectors representing the row, column and data. Alternatively, the function "spconvert" uses a three column matrix format to allow easy importation of data from elsewhere. Created and then filled The function "sparse" or "spalloc" can be used to create an empty matrix that is then filled by the user From a user binary program The user can directly create the sparse matrix within an oct-file. There are several basic functions to return specific sparse matrices. For example the sparse identity matrix, is a matrix that is often needed. It therefore has its own function to create it as 'speye (N)' or 'speye (R, C)', which creates an N-by-N or R-by-C sparse identity matrix. Another typical sparse matrix that is often needed is a random distribution of random elements. The functions "sprand" and "sprandn" perform this for uniform and normal random distributions of elements. They have exactly the same calling convention, where 'sprand (R, C, D)', creates an R-by-C sparse matrix with a density of filled elements of D. Other functions of interest that directly create sparse matrices, are "diag" or its generalization "spdiags", that can take the definition of the diagonals of the matrix and create the sparse matrix that corresponds to this. For example, s = diag (sparse (randn (1,n)), -1); creates a sparse (N+1)-by-(N+1) sparse matrix with a single diagonal defined. -- Function File: B = spdiags (A) -- Function File: [B, D] = spdiags (A) -- Function File: B = spdiags (A, D) -- Function File: A = spdiags (V, D, A) -- Function File: A = spdiags (V, D, M, N) A generalization of the function 'diag'. Called with a single input argument, the nonzero diagonals D of A are extracted. With two arguments the diagonals to extract are given by the vector D. The other two forms of 'spdiags' modify the input matrix by replacing the diagonals. They use the columns of V to replace the diagonals represented by the vector D. If the sparse matrix A is defined then the diagonals of this matrix are replaced. Otherwise a matrix of M by N is created with the diagonals given by the columns of V. Negative values of D represent diagonals below the main diagonal, and positive values of D diagonals above the main diagonal. For example: spdiags (reshape (1:12, 4, 3), [-1 0 1], 5, 4) => 5 10 0 0 1 6 11 0 0 2 7 12 0 0 3 8 0 0 0 4 See also: *note diag: XREFdiag. -- Function File: S = speye (M, N) -- Function File: S = speye (M) -- Function File: S = speye (SZ) Return a sparse identity matrix of size MxN. The implementation is significantly more efficient than 'sparse (eye (M))' as the full matrix is not constructed. Called with a single argument a square matrix of size M-by-M is created. If called with a single vector argument SZ, this argument is taken to be the size of the matrix to create. See also: *note sparse: XREFsparse, *note spdiags: XREFspdiags, *note eye: XREFeye. -- Function File: R = spones (S) Replace the nonzero entries of S with ones. This creates a sparse matrix with the same structure as S. See also: *note sparse: XREFsparse, *note sprand: XREFsprand, *note sprandn: XREFsprandn, *note sprandsym: XREFsprandsym, *note spfun: XREFspfun, *note spy: XREFspy. -- Function File: sprand (M, N, D) -- Function File: sprand (M, N, D, RC) -- Function File: sprand (S) Generate a sparse matrix with uniformly distributed random values. The size of the matrix is MxN with a density of values D. D must be between 0 and 1. Values will be uniformly distributed on the interval (0, 1). If called with a single matrix argument, a sparse matrix is generated with random values wherever the matrix S is nonzero. If called with a scalar fourth argument RC, a random sparse matrix with reciprocal condition number RC is generated. If RC is a vector, then it specifies the first singular values of the generated matrix ('length (RC) <= min (M, N)'). See also: *note sprandn: XREFsprandn, *note sprandsym: XREFsprandsym, *note rand: XREFrand. -- Function File: sprandn (M, N, D) -- Function File: sprandn (M, N, D, RC) -- Function File: sprandn (S) Generate a sparse matrix with normally distributed random values. The size of the matrix is MxN with a density of values D. D must be between 0 and 1. Values will be normally distributed with a mean of 0 and a variance of 1. If called with a single matrix argument, a sparse matrix is generated with random values wherever the matrix S is nonzero. If called with a scalar fourth argument RC, a random sparse matrix with reciprocal condition number RC is generated. If RC is a vector, then it specifies the first singular values of the generated matrix ('length (RC) <= min (M, N)'). See also: *note sprand: XREFsprand, *note sprandsym: XREFsprandsym, *note randn: XREFrandn. -- Function File: sprandsym (N, D) -- Function File: sprandsym (S) Generate a symmetric random sparse matrix. The size of the matrix will be NxN, with a density of values given by D. D must be between 0 and 1 inclusive. Values will be normally distributed with a mean of zero and a variance of 1. If called with a single matrix argument, a random sparse matrix is generated wherever the matrix S is nonzero in its lower triangular part. See also: *note sprand: XREFsprand, *note sprandn: XREFsprandn, *note spones: XREFspones, *note sparse: XREFsparse. The recommended way for the user to create a sparse matrix, is to create two vectors containing the row and column index of the data and a third vector of the same size containing the data to be stored. For example, ri = ci = d = []; for j = 1:c ri = [ri; randperm(r,n)']; ci = [ci; j*ones(n,1)]; d = [d; rand(n,1)]; endfor s = sparse (ri, ci, d, r, c); creates an R-by-C sparse matrix with a random distribution of N ( Compressed Column Sparse (rows=4, cols=4, nnz=3) (1 , 1) -> 1 (2 , 3) -> 2 (3 , 4) -> 3 An example of creating and filling a matrix might be k = 5; nz = r * k; s = spalloc (r, c, nz) for j = 1:c idx = randperm (r); s (:, j) = [zeros(r - k, 1); ... rand(k, 1)] (idx); endfor It should be noted, that due to the way that the Octave assignment functions are written that the assignment will reallocate the memory used by the sparse matrix at each iteration of the above loop. Therefore the "spalloc" function ignores the NZ argument and does not pre-assign the memory for the matrix. Therefore, it is vitally important that code using to above structure should be vectorized as much as possible to minimize the number of assignments and reduce the number of memory allocations. -- Built-in Function: FM = full (SM) Return a full storage matrix from a sparse, diagonal, or permutation matrix, or a range. See also: *note sparse: XREFsparse, *note issparse: XREFissparse. -- Built-in Function: S = spalloc (M, N, NZ) Create an M-by-N sparse matrix with pre-allocated space for at most NZ nonzero elements. This is useful for building a matrix incrementally by a sequence of indexed assignments. Subsequent indexed assignments after 'spalloc' will reuse the pre-allocated memory, provided they are of one of the simple forms * 'S(I:J) = X' * 'S(:,I:J) = X' * 'S(K:L,I:J) = X' and that the following conditions are met: * the assignment does not decrease nnz (S). * after the assignment, nnz (S) does not exceed NZ. * no index is out of bounds. Partial movement of data may still occur, but in general the assignment will be more memory and time efficient under these circumstances. In particular, it is possible to efficiently build a pre-allocated sparse matrix from a contiguous block of columns. The amount of pre-allocated memory for a given matrix may be queried using the function 'nzmax'. See also: *note nzmax: XREFnzmax, *note sparse: XREFsparse. -- Built-in Function: S = sparse (A) -- Built-in Function: S = sparse (I, J, SV, M, N) -- Built-in Function: S = sparse (I, J, SV) -- Built-in Function: S = sparse (M, N) -- Built-in Function: S = sparse (I, J, S, M, N, "unique") -- Built-in Function: S = sparse (I, J, SV, M, N, NZMAX) Create a sparse matrix from a full matrix, or row, column, value triplets. If A is a full matrix, convert it to a sparse matrix representation, removing all zero values in the process. Given the integer index vectors I and J, and a 1-by-'nnz' vector of real or complex values SV, construct the sparse matrix 'S(I(K),J(K)) = SV(K)' with overall dimensions M and N. If any of SV, I or J are scalars, they are expanded to have a common size. If M or N are not specified their values are derived from the maximum index in the vectors I and J as given by 'M = max (I)', 'N = max (J)'. *Note*: if multiple values are specified with the same I, J indices, the corresponding value in S will be the sum of the values at the repeated location. See 'accumarray' for an example of how to produce different behavior, such as taking the minimum instead. If the option "unique" is given, and more than one value is specified at the same I, J indices, then the last specified value will be used. 'sparse (M, N)' will create an empty MxN sparse matrix and is equivalent to 'sparse ([], [], [], M, N)' The argument 'nzmax' is ignored but accepted for compatibility with MATLAB. Example 1 (sum at repeated indices): I = [1 1 2]; J = [1 1 2]; SV = [3 4 5]; sparse (I, J, SV, 3, 4) => Compressed Column Sparse (rows = 3, cols = 4, nnz = 2 [17%]) (1, 1) -> 7 (2, 2) -> 5 Example 2 ("unique" option): I = [1 1 2]; J = [1 1 2]; SV = [3 4 5]; sparse (I, J, SV, 3, 4, "unique") => Compressed Column Sparse (rows = 3, cols = 4, nnz = 2 [17%]) (1, 1) -> 4 (2, 2) -> 5 See also: *note full: XREFfull, *note accumarray: XREFaccumarray, *note spalloc: XREFspalloc, *note spdiags: XREFspdiags, *note speye: XREFspeye, *note spones: XREFspones, *note sprand: XREFsprand, *note sprandn: XREFsprandn, *note sprandsym: XREFsprandsym, *note spconvert: XREFspconvert, *note spfun: XREFspfun. -- Function File: X = spconvert (M) Convert a simple sparse matrix format easily generated by other programs into Octave's internal sparse format. The input M is either a 3 or 4 column real matrix, containing the row, column, real, and imaginary parts of the elements of the sparse matrix. An element with a zero real and imaginary part can be used to force a particular matrix size. See also: *note sparse: XREFsparse. The above problem of memory reallocation can be avoided in oct-files. However, the construction of a sparse matrix from an oct-file is more complex than can be discussed here. *Note External Code Interface::, for a full description of the techniques involved.  File: octave.info, Node: Information, Next: Operators and Functions, Prev: Creating Sparse Matrices, Up: Basics 22.1.3 Finding Information about Sparse Matrices ------------------------------------------------ There are a number of functions that allow information concerning sparse matrices to be obtained. The most basic of these is "issparse" that identifies whether a particular Octave object is in fact a sparse matrix. Another very basic function is "nnz" that returns the number of nonzero entries there are in a sparse matrix, while the function "nzmax" returns the amount of storage allocated to the sparse matrix. Note that Octave tends to crop unused memory at the first opportunity for sparse objects. There are some cases of user created sparse objects where the value returned by "nzmax" will not be the same as "nnz", but in general they will give the same result. The function "spstats" returns some basic statistics on the columns of a sparse matrix including the number of elements, the mean and the variance of each column. -- Built-in Function: issparse (X) Return true if X is a sparse matrix. See also: *note ismatrix: XREFismatrix. -- Built-in Function: N = nnz (A) Return the number of nonzero elements in A. See also: *note nzmax: XREFnzmax, *note nonzeros: XREFnonzeros, *note find: XREFfind. -- Function File: nonzeros (S) Return a vector of the nonzero values of the sparse matrix S. See also: *note find: XREFfind, *note nnz: XREFnnz. -- Built-in Function: N = nzmax (SM) Return the amount of storage allocated to the sparse matrix SM. Note that Octave tends to crop unused memory at the first opportunity for sparse objects. Thus, in general the value of 'nzmax' will be the same as 'nnz' except for some cases of user-created sparse objects. See also: *note nnz: XREFnnz, *note spalloc: XREFspalloc, *note sparse: XREFsparse. -- Function File: [COUNT, MEAN, VAR] = spstats (S) -- Function File: [COUNT, MEAN, VAR] = spstats (S, J) Return the stats for the nonzero elements of the sparse matrix S. COUNT is the number of nonzeros in each column, MEAN is the mean of the nonzeros in each column, and VAR is the variance of the nonzeros in each column. Called with two input arguments, if S is the data and J is the bin number for the data, compute the stats for each bin. In this case, bins can contain data values of zero, whereas with 'spstats (S)' the zeros may disappear. When solving linear equations involving sparse matrices Octave determines the means to solve the equation based on the type of the matrix (*note Sparse Linear Algebra::). Octave probes the matrix type when the div (/) or ldiv (\) operator is first used with the matrix and then caches the type. However the "matrix_type" function can be used to determine the type of the sparse matrix prior to use of the div or ldiv operators. For example, a = tril (sprandn (1024, 1024, 0.02), -1) ... + speye (1024); matrix_type (a); ans = Lower shows that Octave correctly determines the matrix type for lower triangular matrices. "matrix_type" can also be used to force the type of a matrix to be a particular type. For example: a = matrix_type (tril (sprandn (1024, ... 1024, 0.02), -1) + speye (1024), "Lower"); This allows the cost of determining the matrix type to be avoided. However, incorrectly defining the matrix type will result in incorrect results from solutions of linear equations, and so it is entirely the responsibility of the user to correctly identify the matrix type There are several graphical means of finding out information about sparse matrices. The first is the "spy" command, which displays the structure of the nonzero elements of the matrix. *Note Figure 22.1: fig:spmatrix, for an example of the use of "spy". More advanced graphical information can be obtained with the "treeplot", "etreeplot" and "gplot" commands. [image src="spmatrix.png" text=" | * * | * * * * | * * * * | * * * * 5 - * * * * | * * * * | * * * * | * * * | * * 10 - * * | * * | * * | * * | * * 15 - * * |----------|---------|---------| 5 10 15"] Figure 22.1: Structure of simple sparse matrix. One use of sparse matrices is in graph theory, where the interconnections between nodes are represented as an adjacency matrix. That is, if the i-th node in a graph is connected to the j-th node. Then the ij-th node (and in the case of undirected graphs the ji-th node) of the sparse adjacency matrix is nonzero. If each node is then associated with a set of coordinates, then the "gplot" command can be used to graphically display the interconnections between nodes. As a trivial example of the use of "gplot" consider the example, A = sparse ([2,6,1,3,2,4,3,5,4,6,1,5], [1,1,2,2,3,3,4,4,5,5,6,6],1,6,6); xy = [0,4,8,6,4,2;5,0,5,7,5,7]'; gplot (A,xy) which creates an adjacency matrix 'A' where node 1 is connected to nodes 2 and 6, node 2 with nodes 1 and 3, etc. The coordinates of the nodes are given in the n-by-2 matrix 'xy'. The dependencies between the nodes of a Cholesky factorization can be calculated in linear time without explicitly needing to calculate the Cholesky factorization by the 'etree' command. This command returns the elimination tree of the matrix and can be displayed graphically by the command 'treeplot (etree (A))' if 'A' is symmetric or 'treeplot (etree (A+A'))' otherwise. -- Function File: spy (X) -- Function File: spy (..., MARKERSIZE) -- Function File: spy (..., LINE_SPEC) Plot the sparsity pattern of the sparse matrix X. If the argument MARKERSIZE is given as a scalar value, it is used to determine the point size in the plot. If the string LINE_SPEC is given it is passed to 'plot' and determines the appearance of the plot. See also: *note plot: XREFplot, *note gplot: XREFgplot. -- Loadable Function: P = etree (S) -- Loadable Function: P = etree (S, TYP) -- Loadable Function: [P, Q] = etree (S, TYP) Return the elimination tree for the matrix S. By default S is assumed to be symmetric and the symmetric elimination tree is returned. The argument TYP controls whether a symmetric or column elimination tree is returned. Valid values of TYP are "sym" or "col", for symmetric or column elimination tree respectively. Called with a second argument, 'etree' also returns the postorder permutations on the tree. -- Function File: etreeplot (A) -- Function File: etreeplot (A, NODE_STYLE, EDGE_STYLE) Plot the elimination tree of the matrix A or A+A' if A in not symmetric. The optional parameters NODE_STYLE and EDGE_STYLE define the output style. See also: *note treeplot: XREFtreeplot, *note gplot: XREFgplot. -- Function File: gplot (A, XY) -- Function File: gplot (A, XY, LINE_STYLE) -- Function File: [X, Y] = gplot (A, XY) Plot a graph defined by A and XY in the graph theory sense. A is the adjacency matrix of the array to be plotted and XY is an N-by-2 matrix containing the coordinates of the nodes of the graph. The optional parameter LINE_STYLE defines the output style for the plot. Called with no output arguments the graph is plotted directly. Otherwise, return the coordinates of the plot in X and Y. See also: *note treeplot: XREFtreeplot, *note etreeplot: XREFetreeplot, *note spy: XREFspy. -- Function File: treeplot (TREE) -- Function File: treeplot (TREE, NODE_STYLE, EDGE_STYLE) Produce a graph of tree or forest. The first argument is vector of predecessors. The optional parameters NODE_STYLE and EDGE_STYLE define the output plot style. The complexity of the algorithm is O(n) in terms of is time and memory requirements. See also: *note etreeplot: XREFetreeplot, *note gplot: XREFgplot. -- Function File: treelayout (TREE) -- Function File: treelayout (TREE, PERMUTATION) treelayout lays out a tree or a forest. The first argument TREE is a vector of predecessors. The parameter PERMUTATION is an optional postorder permutation. The complexity of the algorithm is O(n) in terms of time and memory requirements. See also: *note etreeplot: XREFetreeplot, *note gplot: XREFgplot, *note treeplot: XREFtreeplot.  File: octave.info, Node: Operators and Functions, Prev: Information, Up: Basics 22.1.4 Basic Operators and Functions on Sparse Matrices ------------------------------------------------------- * Menu: * Sparse Functions:: * Return Types of Operators and Functions:: * Mathematical Considerations::  File: octave.info, Node: Sparse Functions, Next: Return Types of Operators and Functions, Up: Operators and Functions 22.1.4.1 Sparse Functions ......................... Many Octave functions have been overloaded to work with either sparse or full matrices. There is no difference in calling convention when using an overloaded function with a sparse matrix, however, there is also no access to potentially sparse-specific features. At any time the sparse matrix specific version of a function can be used by explicitly calling its function name. The table below lists all of the sparse functions of Octave. Note that the names of the specific sparse forms of the functions are typically the same as the general versions with a "sp" prefix. In the table below, and in the rest of this article, the specific sparse versions of functions are used. Generate sparse matrices: "spalloc", "spdiags", "speye", "sprand", "sprandn", "sprandsym" Sparse matrix conversion: "full", "sparse", "spconvert" Manipulate sparse matrices "issparse", "nnz", "nonzeros", "nzmax", "spfun", "spones", "spy" Graph Theory: "etree", "etreeplot", "gplot", "treeplot" Sparse matrix reordering: "amd", "ccolamd", "colamd", "colperm", "csymamd", "dmperm", "symamd", "randperm", "symrcm" Linear algebra: "condest", "eigs", "matrix_type", "normest", "sprank", "spaugment", "svds" Iterative techniques: "ichol", "ilu", "pcg", "pcr" Miscellaneous: "spparms", "symbfact", "spstats" In addition all of the standard Octave mapper functions (i.e., basic math functions that take a single argument) such as "abs", etc. can accept sparse matrices. The reader is referred to the documentation supplied with these functions within Octave itself for further details.  File: octave.info, Node: Return Types of Operators and Functions, Next: Mathematical Considerations, Prev: Sparse Functions, Up: Operators and Functions 22.1.4.2 Return Types of Operators and Functions ................................................ The two basic reasons to use sparse matrices are to reduce the memory usage and to not have to do calculations on zero elements. The two are closely related in that the computation time on a sparse matrix operator or function is roughly linear with the number of nonzero elements. Therefore, there is a certain density of nonzero elements of a matrix where it no longer makes sense to store it as a sparse matrix, but rather as a full matrix. For this reason operators and functions that have a high probability of returning a full matrix will always return one. For example adding a scalar constant to a sparse matrix will almost always make it a full matrix, and so the example, speye (3) + 0 => 1 0 0 0 1 0 0 0 1 returns a full matrix as can be seen. Additionally, if 'sparse_auto_mutate' is true, all sparse functions test the amount of memory occupied by the sparse matrix to see if the amount of storage used is larger than the amount used by the full equivalent. Therefore 'speye (2) * 1' will return a full matrix as the memory used is smaller for the full version than the sparse version. As all of the mixed operators and functions between full and sparse matrices exist, in general this does not cause any problems. However, one area where it does cause a problem is where a sparse matrix is promoted to a full matrix, where subsequent operations would resparsify the matrix. Such cases are rare, but can be artificially created, for example '(fliplr (speye (3)) + speye (3)) - speye (3)' gives a full matrix when it should give a sparse one. In general, where such cases occur, they impose only a small memory penalty. There is however one known case where this behavior of Octave's sparse matrices will cause a problem. That is in the handling of the "diag" function. Whether "diag" returns a sparse or full matrix depending on the type of its input arguments. So a = diag (sparse ([1,2,3]), -1); should return a sparse matrix. To ensure this actually happens, the "sparse" function, and other functions based on it like "speye", always returns a sparse matrix, even if the memory used will be larger than its full representation. -- Built-in Function: VAL = sparse_auto_mutate () -- Built-in Function: OLD_VAL = sparse_auto_mutate (NEW_VAL) -- Built-in Function: sparse_auto_mutate (NEW_VAL, "local") Query or set the internal variable that controls whether Octave will automatically mutate sparse matrices to full matrices to save memory. For example: s = speye (3); sparse_auto_mutate (false); s(:, 1) = 1; typeinfo (s) => sparse matrix sparse_auto_mutate (true); s(1, :) = 1; typeinfo (s) => matrix When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. Note that the 'sparse_auto_mutate' option is incompatible with MATLAB, and so it is off by default.  File: octave.info, Node: Mathematical Considerations, Prev: Return Types of Operators and Functions, Up: Operators and Functions 22.1.4.3 Mathematical Considerations .................................... The attempt has been made to make sparse matrices behave in exactly the same manner as there full counterparts. However, there are certain differences and especially differences with other products sparse implementations. First, the "./" and ".^" operators must be used with care. Consider what the examples s = speye (4); a1 = s .^ 2; a2 = s .^ s; a3 = s .^ -2; a4 = s ./ 2; a5 = 2 ./ s; a6 = s ./ s; will give. The first example of S raised to the power of 2 causes no problems. However S raised element-wise to itself involves a large number of terms '0 .^ 0' which is 1. There 'S .^ S' is a full matrix. Likewise 'S .^ -2' involves terms like '0 .^ -2' which is infinity, and so 'S .^ -2' is equally a full matrix. For the "./" operator 'S ./ 2' has no problems, but '2 ./ S' involves a large number of infinity terms as well and is equally a full matrix. The case of 'S ./ S' involves terms like '0 ./ 0' which is a 'NaN' and so this is equally a full matrix with the zero elements of S filled with 'NaN' values. The above behavior is consistent with full matrices, but is not consistent with sparse implementations in other products. A particular problem of sparse matrices comes about due to the fact that as the zeros are not stored, the sign-bit of these zeros is equally not stored. In certain cases the sign-bit of zero is important. For example: a = 0 ./ [-1, 1; 1, -1]; b = 1 ./ a => -Inf Inf Inf -Inf c = 1 ./ sparse (a) => Inf Inf Inf Inf To correct this behavior would mean that zero elements with a negative sign-bit would need to be stored in the matrix to ensure that their sign-bit was respected. This is not done at this time, for reasons of efficiency, and so the user is warned that calculations where the sign-bit of zero is important must not be done using sparse matrices. In general any function or operator used on a sparse matrix will result in a sparse matrix with the same or a larger number of nonzero elements than the original matrix. This is particularly true for the important case of sparse matrix factorizations. The usual way to address this is to reorder the matrix, such that its factorization is sparser than the factorization of the original matrix. That is the factorization of 'L * U = P * S * Q' has sparser terms 'L' and 'U' than the equivalent factorization 'L * U = S'. Several functions are available to reorder depending on the type of the matrix to be factorized. If the matrix is symmetric positive-definite, then "symamd" or "csymamd" should be used. Otherwise "amd", "colamd" or "ccolamd" should be used. For completeness the reordering functions "colperm" and "randperm" are also available. *Note Figure 22.2: fig:simplematrix, for an example of the structure of a simple positive definite matrix. [image src="spmatrix.png" text=" | * * | * * * * | * * * * | * * * * 5 - * * * * | * * * * | * * * * | * * * | * * 10 - * * | * * | * * | * * | * * 15 - * * |----------|---------|---------| 5 10 15"] Figure 22.2: Structure of simple sparse matrix. The standard Cholesky factorization of this matrix can be obtained by the same command that would be used for a full matrix. This can be visualized with the command 'r = chol (A); spy (r);'. *Note Figure 22.3: fig:simplechol. The original matrix had 43 nonzero terms, while this Cholesky factorization has 71, with only half of the symmetric matrix being stored. This is a significant level of fill in, and although not an issue for such a small test case, can represents a large overhead in working with other sparse matrices. The appropriate sparsity preserving permutation of the original matrix is given by "symamd" and the factorization using this reordering can be visualized using the command 'q = symamd (A); r = chol (A(q,q)); spy (r)'. This gives 29 nonzero terms which is a significant improvement. The Cholesky factorization itself can be used to determine the appropriate sparsity preserving reordering of the matrix during the factorization, In that case this might be obtained with three return arguments as '[r, p, q] = chol (A); spy (r)'. [image src="spchol.png" text=" | * * | * * * | * * * * | * * * * * 5 - * * * * * * | * * * * * * * | * * * * * * * * | * * * * * * * * | * * * * * * * 10 - * * * * * * | * * * * * | * * * * | * * * | * * 15 - * |----------|---------|---------| 5 10 15"] Figure 22.3: Structure of the unpermuted Cholesky factorization of the above matrix. [image src="spcholperm.png" text=" | * * | * * | * * | * * 5 - * * | * * | * * | * * | * * 10 - * * | * * | * * | * * | * * 15 - * |----------|---------|---------| 5 10 15"] Figure 22.4: Structure of the permuted Cholesky factorization of the above matrix. In the case of an asymmetric matrix, the appropriate sparsity preserving permutation is "colamd" and the factorization using this reordering can be visualized using the command 'q = colamd (A); [l, u, p] = lu (A(:,q)); spy (l+u)'. Finally, Octave implicitly reorders the matrix when using the div (/) and ldiv (\) operators, and so no the user does not need to explicitly reorder the matrix to maximize performance. -- Loadable Function: P = amd (S) -- Loadable Function: P = amd (S, OPTS) Return the approximate minimum degree permutation of a matrix. This is a permutation such that the Cholesky factorization of 'S (P, P)' tends to be sparser than the Cholesky factorization of S itself. 'amd' is typically faster than 'symamd' but serves a similar purpose. The optional parameter OPTS is a structure that controls the behavior of 'amd'. The fields of the structure are OPTS.dense Determines what 'amd' considers to be a dense row or column of the input matrix. Rows or columns with more than 'max (16, (dense * sqrt (N)))' entries, where N is the order of the matrix S, are ignored by 'amd' during the calculation of the permutation. The value of dense must be a positive scalar and the default value is 10.0 OPTS.aggressive If this value is a nonzero scalar, then 'amd' performs aggressive absorption. The default is not to perform aggressive absorption. The author of the code itself is Timothy A. Davis , University of Florida (see ). See also: *note symamd: XREFsymamd, *note colamd: XREFcolamd. -- Loadable Function: P = ccolamd (S) -- Loadable Function: P = ccolamd (S, KNOBS) -- Loadable Function: P = ccolamd (S, KNOBS, CMEMBER) -- Loadable Function: [P, STATS] = ccolamd (...) Constrained column approximate minimum degree permutation. 'P = ccolamd (S)' returns the column approximate minimum degree permutation vector for the sparse matrix S. For a non-symmetric matrix S, 'S(:, P)' tends to have sparser LU factors than S. 'chol (S(:, P)' * S(:, P))' also tends to be sparser than 'chol (S' * S)'. 'P = ccolamd (S, 1)' optimizes the ordering for 'lu (S(:, P))'. The ordering is followed by a column elimination tree post-ordering. KNOBS is an optional 1-element to 5-element input vector, with a default value of '[0 10 10 1 0]' if not present or empty. Entries not present are set to their defaults. 'KNOBS(1)' if nonzero, the ordering is optimized for 'lu (S(:, p))'. It will be a poor ordering for 'chol (S(:, P)' * S(:, P))'. This is the most important knob for ccolamd. 'KNOBS(2)' if S is m-by-n, rows with more than 'max (16, KNOBS(2) * sqrt (n))' entries are ignored. 'KNOBS(3)' columns with more than 'max (16, KNOBS(3) * sqrt (min (M, N)))' entries are ignored and ordered last in the output permutation (subject to the cmember constraints). 'KNOBS(4)' if nonzero, aggressive absorption is performed. 'KNOBS(5)' if nonzero, statistics and knobs are printed. CMEMBER is an optional vector of length n. It defines the constraints on the column ordering. If 'CMEMBER(j) = C', then column J is in constraint set C (C must be in the range 1 to n). In the output permutation P, all columns in set 1 appear first, followed by all columns in set 2, and so on. 'CMEMBER = ones (1,n)' if not present or empty. 'ccolamd (S, [], 1 : n)' returns '1 : n' 'P = ccolamd (S)' is about the same as 'P = colamd (S)'. KNOBS and its default values differ. 'colamd' always does aggressive absorption, and it finds an ordering suitable for both 'lu (S(:, P))' and 'chol (S(:, P)' * S(:, P))'; it cannot optimize its ordering for 'lu (S(:, P))' to the extent that 'ccolamd (S, 1)' can. STATS is an optional 20-element output vector that provides data about the ordering and the validity of the input matrix S. Ordering statistics are in 'STATS(1 : 3)'. 'STATS(1)' and 'STATS(2)' are the number of dense or empty rows and columns ignored by CCOLAMD and 'STATS(3)' is the number of garbage collections performed on the internal data structure used by CCOLAMD (roughly of size '2.2 * nnz (S) + 4 * M + 7 * N' integers). 'STATS(4 : 7)' provide information if CCOLAMD was able to continue. The matrix is OK if 'STATS(4)' is zero, or 1 if invalid. 'STATS(5)' is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. 'STATS(6)' is the last seen duplicate or out-of-order row index in the column index given by 'STATS(5)', or zero if no such row index exists. 'STATS(7)' is the number of duplicate or out-of-order row indices. 'STATS(8 : 20)' is always zero in the current version of CCOLAMD (reserved for future use). The authors of the code itself are S. Larimore, T. Davis (Univ. of Florida) and S. Rajamanickam in collaboration with J. Bilbert and E. Ng. Supported by the National Science Foundation (DMS-9504974, DMS-9803599, CCR-0203270), and a grant from Sandia National Lab. See for ccolamd, csymamd, amd, colamd, symamd, and other related orderings. See also: *note colamd: XREFcolamd, *note csymamd: XREFcsymamd. -- Loadable Function: P = colamd (S) -- Loadable Function: P = colamd (S, KNOBS) -- Loadable Function: [P, STATS] = colamd (S) -- Loadable Function: [P, STATS] = colamd (S, KNOBS) Compute the column approximate minimum degree permutation. 'P = colamd (S)' returns the column approximate minimum degree permutation vector for the sparse matrix S. For a non-symmetric matrix S, 'S(:,P)' tends to have sparser LU factors than S. The Cholesky factorization of 'S(:,P)' * S(:,P)' also tends to be sparser than that of 'S' * S'. KNOBS is an optional one- to three-element input vector. If S is m-by-n, then rows with more than 'max(16,KNOBS(1)*sqrt(n))' entries are ignored. Columns with more than 'max (16,KNOBS(2)*sqrt(min(m,n)))' entries are removed prior to ordering, and ordered last in the output permutation P. Only completely dense rows or columns are removed if 'KNOBS(1)' and 'KNOBS(2)' are < 0, respectively. If 'KNOBS(3)' is nonzero, STATS and KNOBS are printed. The default is 'KNOBS = [10 10 0]'. Note that KNOBS differs from earlier versions of colamd. STATS is an optional 20-element output vector that provides data about the ordering and the validity of the input matrix S. Ordering statistics are in 'STATS(1:3)'. 'STATS(1)' and 'STATS(2)' are the number of dense or empty rows and columns ignored by COLAMD and 'STATS(3)' is the number of garbage collections performed on the internal data structure used by COLAMD (roughly of size '2.2 * nnz(S) + 4 * M + 7 * N' integers). Octave built-in functions are intended to generate valid sparse matrices, with no duplicate entries, with ascending row indices of the nonzeros in each column, with a non-negative number of entries in each column (!) and so on. If a matrix is invalid, then COLAMD may or may not be able to continue. If there are duplicate entries (a row index appears two or more times in the same column) or if the row indices in a column are out of order, then COLAMD can correct these errors by ignoring the duplicate entries and sorting each column of its internal copy of the matrix S (the input matrix S is not repaired, however). If a matrix is invalid in other ways then COLAMD cannot continue, an error message is printed, and no output arguments (P or STATS) are returned. COLAMD is thus a simple way to check a sparse matrix to see if it's valid. 'STATS(4:7)' provide information if COLAMD was able to continue. The matrix is OK if 'STATS(4)' is zero, or 1 if invalid. 'STATS(5)' is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. 'STATS(6)' is the last seen duplicate or out-of-order row index in the column index given by 'STATS(5)', or zero if no such row index exists. 'STATS(7)' is the number of duplicate or out-of-order row indices. 'STATS(8:20)' is always zero in the current version of COLAMD (reserved for future use). The ordering is followed by a column elimination tree post-ordering. The authors of the code itself are Stefan I. Larimore and Timothy A. Davis , University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. (see ) See also: *note colperm: XREFcolperm, *note symamd: XREFsymamd, *note ccolamd: XREFccolamd. -- Function File: P = colperm (S) Return the column permutations such that the columns of 'S (:, P)' are ordered in terms of increasing number of nonzero elements. If S is symmetric, then P is chosen such that 'S (P, P)' orders the rows and columns with increasing number of nonzeros elements. -- Loadable Function: P = csymamd (S) -- Loadable Function: P = csymamd (S, KNOBS) -- Loadable Function: P = csymamd (S, KNOBS, CMEMBER) -- Loadable Function: [P, STATS] = csymamd (...) For a symmetric positive definite matrix S, return the permutation vector P such that 'S(P,P)' tends to have a sparser Cholesky factor than S. Sometimes 'csymamd' works well for symmetric indefinite matrices too. The matrix S is assumed to be symmetric; only the strictly lower triangular part is referenced. S must be square. The ordering is followed by an elimination tree post-ordering. KNOBS is an optional 1-element to 3-element input vector, with a default value of '[10 1 0]'. Entries not present are set to their defaults. 'KNOBS(1)' If S is n-by-n, then rows and columns with more than 'max(16,KNOBS(1)*sqrt(n))' entries are ignored, and ordered last in the output permutation (subject to the cmember constraints). 'KNOBS(2)' If nonzero, aggressive absorption is performed. 'KNOBS(3)' If nonzero, statistics and knobs are printed. CMEMBER is an optional vector of length n. It defines the constraints on the ordering. If 'CMEMBER(j) = S', then row/column j is in constraint set C (C must be in the range 1 to n). In the output permutation P, rows/columns in set 1 appear first, followed by all rows/columns in set 2, and so on. 'CMEMBER = ones (1,n)' if not present or empty. 'csymamd (S,[],1:n)' returns '1:n'. 'P = csymamd (S)' is about the same as 'P = symamd (S)'. KNOBS and its default values differ. 'STATS(4:7)' provide information if CCOLAMD was able to continue. The matrix is OK if 'STATS(4)' is zero, or 1 if invalid. 'STATS(5)' is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. 'STATS(6)' is the last seen duplicate or out-of-order row index in the column index given by 'STATS(5)', or zero if no such row index exists. 'STATS(7)' is the number of duplicate or out-of-order row indices. 'STATS(8:20)' is always zero in the current version of CCOLAMD (reserved for future use). The authors of the code itself are S. Larimore, T. Davis (Univ. of Florida) and S. Rajamanickam in collaboration with J. Bilbert and E. Ng. Supported by the National Science Foundation (DMS-9504974, DMS-9803599, CCR-0203270), and a grant from Sandia National Lab. See for ccolamd, csymamd, amd, colamd, symamd, and other related orderings. See also: *note symamd: XREFsymamd, *note ccolamd: XREFccolamd. -- Loadable Function: P = dmperm (S) -- Loadable Function: [P, Q, R, S] = dmperm (S) Perform a Dulmage-Mendelsohn permutation of the sparse matrix S. With a single output argument 'dmperm' performs the row permutations P such that 'S(P,:)' has no zero elements on the diagonal. Called with two or more output arguments, returns the row and column permutations, such that 'S(P, Q)' is in block triangular form. The values of R and S define the boundaries of the blocks. If S is square then 'R == S'. The method used is described in: A. Pothen & C.-J. Fan. 'Computing the Block Triangular Form of a Sparse Matrix'. ACM Trans. Math. Software, 16(4):303-324, 1990. See also: *note colamd: XREFcolamd, *note ccolamd: XREFccolamd. -- Loadable Function: P = symamd (S) -- Loadable Function: P = symamd (S, KNOBS) -- Loadable Function: [P, STATS] = symamd (S) -- Loadable Function: [P, STATS] = symamd (S, KNOBS) For a symmetric positive definite matrix S, returns the permutation vector p such that 'S(P, P)' tends to have a sparser Cholesky factor than S. Sometimes 'symamd' works well for symmetric indefinite matrices too. The matrix S is assumed to be symmetric; only the strictly lower triangular part is referenced. S must be square. KNOBS is an optional one- to two-element input vector. If S is n-by-n, then rows and columns with more than 'max (16,KNOBS(1)*sqrt(n))' entries are removed prior to ordering, and ordered last in the output permutation P. No rows/columns are removed if 'KNOBS(1) < 0'. If 'KNOBS (2)' is nonzero, 'stats' and KNOBS are printed. The default is 'KNOBS = [10 0]'. Note that KNOBS differs from earlier versions of 'symamd'. STATS is an optional 20-element output vector that provides data about the ordering and the validity of the input matrix S. Ordering statistics are in 'STATS(1:3)'. 'STATS(1) = STATS(2)' is the number of dense or empty rows and columns ignored by SYMAMD and 'STATS(3)' is the number of garbage collections performed on the internal data structure used by SYMAMD (roughly of size '8.4 * nnz (tril (S, -1)) + 9 * N' integers). Octave built-in functions are intended to generate valid sparse matrices, with no duplicate entries, with ascending row indices of the nonzeros in each column, with a non-negative number of entries in each column (!) and so on. If a matrix is invalid, then SYMAMD may or may not be able to continue. If there are duplicate entries (a row index appears two or more times in the same column) or if the row indices in a column are out of order, then SYMAMD can correct these errors by ignoring the duplicate entries and sorting each column of its internal copy of the matrix S (the input matrix S is not repaired, however). If a matrix is invalid in other ways then SYMAMD cannot continue, an error message is printed, and no output arguments (P or STATS) are returned. SYMAMD is thus a simple way to check a sparse matrix to see if it's valid. 'STATS(4:7)' provide information if SYMAMD was able to continue. The matrix is OK if 'STATS (4)' is zero, or 1 if invalid. 'STATS(5)' is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. 'STATS(6)' is the last seen duplicate or out-of-order row index in the column index given by 'STATS(5)', or zero if no such row index exists. 'STATS(7)' is the number of duplicate or out-of-order row indices. 'STATS(8:20)' is always zero in the current version of SYMAMD (reserved for future use). The ordering is followed by a column elimination tree post-ordering. The authors of the code itself are Stefan I. Larimore and Timothy A. Davis , University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. (see ) See also: *note colperm: XREFcolperm, *note colamd: XREFcolamd. -- Loadable Function: P = symrcm (S) Return the symmetric reverse Cuthill-McKee permutation of S. P is a permutation vector such that 'S(P, P)' tends to have its diagonal elements closer to the diagonal than S. This is a good preordering for LU or Cholesky factorization of matrices that come from "long, skinny" problems. It works for both symmetric and asymmetric S. The algorithm represents a heuristic approach to the NP-complete bandwidth minimization problem. The implementation is based in the descriptions found in E. Cuthill, J. McKee. 'Reducing the Bandwidth of Sparse Symmetric Matrices'. Proceedings of the 24th ACM National Conference, 157-172 1969, Brandon Press, New Jersey. A. George, J.W.H. Liu. 'Computer Solution of Large Sparse Positive Definite Systems', Prentice Hall Series in Computational Mathematics, ISBN 0-13-165274-5, 1981. See also: *note colperm: XREFcolperm, *note colamd: XREFcolamd, *note symamd: XREFsymamd.  File: octave.info, Node: Sparse Linear Algebra, Next: Iterative Techniques, Prev: Basics, Up: Sparse Matrices 22.2 Linear Algebra on Sparse Matrices ====================================== Octave includes a polymorphic solver for sparse matrices, where the exact solver used to factorize the matrix, depends on the properties of the sparse matrix itself. Generally, the cost of determining the matrix type is small relative to the cost of factorizing the matrix itself, but in any case the matrix type is cached once it is calculated, so that it is not re-determined each time it is used in a linear equation. The selection tree for how the linear equation is solve is 1. If the matrix is diagonal, solve directly and goto 8 2. If the matrix is a permuted diagonal, solve directly taking into account the permutations. Goto 8 3. If the matrix is square, banded and if the band density is less than that given by 'spparms ("bandden")' continue, else goto 4. a. If the matrix is tridiagonal and the right-hand side is not sparse continue, else goto 3b. 1. If the matrix is Hermitian, with a positive real diagonal, attempt Cholesky factorization using LAPACK xPTSV. 2. If the above failed or the matrix is not Hermitian with a positive real diagonal use Gaussian elimination with pivoting using LAPACK xGTSV, and goto 8. b. If the matrix is Hermitian with a positive real diagonal, attempt Cholesky factorization using LAPACK xPBTRF. c. if the above failed or the matrix is not Hermitian with a positive real diagonal use Gaussian elimination with pivoting using LAPACK xGBTRF, and goto 8. 4. If the matrix is upper or lower triangular perform a sparse forward or backward substitution, and goto 8 5. If the matrix is an upper triangular matrix with column permutations or lower triangular matrix with row permutations, perform a sparse forward or backward substitution, and goto 8 6. If the matrix is square, Hermitian with a real positive diagonal, attempt sparse Cholesky factorization using CHOLMOD. 7. If the sparse Cholesky factorization failed or the matrix is not Hermitian with a real positive diagonal, and the matrix is square, factorize using UMFPACK. 8. If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a minimum norm solution using CXSPARSE(1). The band density is defined as the number of nonzero values in the band divided by the total number of values in the full band. The banded matrix solvers can be entirely disabled by using "spparms" to set 'bandden' to 1 (i.e., 'spparms ("bandden", 1)'). The QR solver factorizes the problem with a Dulmage-Mendelsohn decomposition, to separate the problem into blocks that can be treated as over-determined, multiple well determined blocks, and a final over-determined block. For matrices with blocks of strongly connected nodes this is a big win as LU decomposition can be used for many blocks. It also significantly improves the chance of finding a solution to over-determined problems rather than just returning a vector of "NaN"'s. All of the solvers above, can calculate an estimate of the condition number. This can be used to detect numerical stability problems in the solution and force a minimum norm solution to be used. However, for narrow banded, triangular or diagonal matrices, the cost of calculating the condition number is significant, and can in fact exceed the cost of factoring the matrix. Therefore the condition number is not calculated in these cases, and Octave relies on simpler techniques to detect singular matrices or the underlying LAPACK code in the case of banded matrices. The user can force the type of the matrix with the 'matrix_type' function. This overcomes the cost of discovering the type of the matrix. However, it should be noted that identifying the type of the matrix incorrectly will lead to unpredictable results, and so 'matrix_type' should be used with care. -- Function File: N = normest (A) -- Function File: N = normest (A, TOL) -- Function File: [N, C] = normest (...) Estimate the 2-norm of the matrix A using a power series analysis. This is typically used for large matrices, where the cost of calculating 'norm (A)' is prohibitive and an approximation to the 2-norm is acceptable. TOL is the tolerance to which the 2-norm is calculated. By default TOL is 1e-6. The optional output C returns the number of iterations needed for 'normest' to converge. -- Function File: [EST, V, W, ITER] = onenormest (A, T) -- Function File: [EST, V, W, ITER] = onenormest (APPLY, APPLY_T, N, T) Apply Higham and Tisseur's randomized block 1-norm estimator to matrix A using T test vectors. If T exceeds 5, then only 5 test vectors are used. If the matrix is not explicit, e.g., when estimating the norm of 'inv (A)' given an LU factorization, 'onenormest' applies A and its conjugate transpose through a pair of functions APPLY and APPLY_T, respectively, to a dense matrix of size N by T. The implicit version requires an explicit dimension N. Returns the norm estimate EST, two vectors V and W related by norm '(W, 1) = EST * norm (V, 1)', and the number of iterations ITER. The number of iterations is limited to 10 and is at least 2. References: * N.J. Higham and F. Tisseur, 'A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra'. SIMAX vol 21, no 4, pp 1185-1201. * N.J. Higham and F. Tisseur, 'A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra'. See also: *note condest: XREFcondest, *note norm: XREFnorm, *note cond: XREFcond. -- Function File: condest (A) -- Function File: condest (A, T) -- Function File: [EST, V] = condest (...) -- Function File: [EST, V] = condest (A, SOLVE, SOLVE_T, T) -- Function File: [EST, V] = condest (APPLY, APPLY_T, SOLVE, SOLVE_T, N, T) Estimate the 1-norm condition number of a matrix A using T test vectors using a randomized 1-norm estimator. If T exceeds 5, then only 5 test vectors are used. If the matrix is not explicit, e.g., when estimating the condition number of A given an LU factorization, 'condest' uses the following functions: APPLY 'A*x' for a matrix 'x' of size N by T. APPLY_T 'A'*x' for a matrix 'x' of size N by T. SOLVE 'A \ b' for a matrix 'b' of size N by T. SOLVE_T 'A' \ b' for a matrix 'b' of size N by T. The implicit version requires an explicit dimension N. 'condest' uses a randomized algorithm to approximate the 1-norms. 'condest' returns the 1-norm condition estimate EST and a vector V satisfying 'norm (A*v, 1) == norm (A, 1) * norm (V, 1) / EST'. When EST is large, V is an approximate null vector. References: * N.J. Higham and F. Tisseur, 'A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra'. SIMAX vol 21, no 4, pp 1185-1201. * N.J. Higham and F. Tisseur, 'A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra'. See also: *note cond: XREFcond, *note norm: XREFnorm, *note onenormest: XREFonenormest. -- Built-in Function: spparms () -- Built-in Function: VALS = spparms () -- Built-in Function: [KEYS, VALS] = spparms () -- Built-in Function: VAL = spparms (KEY) -- Built-in Function: spparms (VALS) -- Built-in Function: spparms ("default") -- Built-in Function: spparms ("tight") -- Built-in Function: spparms (KEY, VAL) Query or set the parameters used by the sparse solvers and factorization functions. The first four calls above get information about the current settings, while the others change the current settings. The parameters are stored as pairs of keys and values, where the values are all floats and the keys are one of the following strings: 'spumoni' Printing level of debugging information of the solvers (default 0) 'ths_rel' Included for compatibility. Not used. (default 1) 'ths_abs' Included for compatibility. Not used. (default 1) 'exact_d' Included for compatibility. Not used. (default 0) 'supernd' Included for compatibility. Not used. (default 3) 'rreduce' Included for compatibility. Not used. (default 3) 'wh_frac' Included for compatibility. Not used. (default 0.5) 'autommd' Flag whether the LU/QR and the '\' and '/' operators will automatically use the sparsity preserving mmd functions (default 1) 'autoamd' Flag whether the LU and the '\' and '/' operators will automatically use the sparsity preserving amd functions (default 1) 'piv_tol' The pivot tolerance of the UMFPACK solvers (default 0.1) 'sym_tol' The pivot tolerance of the UMFPACK symmetric solvers (default 0.001) 'bandden' The density of nonzero elements in a banded matrix before it is treated by the LAPACK banded solvers (default 0.5) 'umfpack' Flag whether the UMFPACK or mmd solvers are used for the LU, '\' and '/' operations (default 1) The value of individual keys can be set with 'spparms (KEY, VAL)'. The default values can be restored with the special keyword "default". The special keyword "tight" can be used to set the mmd solvers to attempt a sparser solution at the potential cost of longer running time. See also: *note chol: XREFchol, *note colamd: XREFcolamd, *note lu: XREFlu, *note qr: XREFqr, *note symamd: XREFsymamd. -- Loadable Function: P = sprank (S) Calculate the structural rank of the sparse matrix S. Note that only the structure of the matrix is used in this calculation based on a Dulmage-Mendelsohn permutation to block triangular form. As such the numerical rank of the matrix S is bounded by 'sprank (S) >= rank (S)'. Ignoring floating point errors 'sprank (S) == rank (S)'. See also: *note dmperm: XREFdmperm. -- Loadable Function: [COUNT, H, PARENT, POST, R] = symbfact (S) -- Loadable Function: [...] = symbfact (S, TYP) -- Loadable Function: [...] = symbfact (S, TYP, MODE) Perform a symbolic factorization analysis on the sparse matrix S. The input variables are S S is a complex or real sparse matrix. TYP Is the type of the factorization and can be one of 'sym' Factorize S. This is the default. 'col' Factorize 'S' * S'. 'row' Factorize S * S'. 'lo' Factorize S' MODE The default is to return the Cholesky factorization for R, and if MODE is 'L', the conjugate transpose of the Cholesky factorization is returned. The conjugate transpose version is faster and uses less memory, but returns the same values for COUNT, H, PARENT and POST outputs. The output variables are COUNT The row counts of the Cholesky factorization as determined by TYP. H The height of the elimination tree. PARENT The elimination tree itself. POST A sparse boolean matrix whose structure is that of the Cholesky factorization as determined by TYP. For non square matrices, the user can also utilize the 'spaugment' function to find a least squares solution to a linear equation. -- Function File: S = spaugment (A, C) Create the augmented matrix of A. This is given by [C * eye(M, M), A; A', zeros(N, N)] This is related to the least squares solution of 'A \ B', by S * [ R / C; x] = [ B, zeros(N, columns(B)) ] where R is the residual error R = B - A * X As the matrix S is symmetric indefinite it can be factorized with 'lu', and the minimum norm solution can therefore be found without the need for a 'qr' factorization. As the residual error will be 'zeros (M, M)' for underdetermined problems, and example can be m = 11; n = 10; mn = max (m, n); A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)], [-1, 0, 1], m, n); x0 = A \ ones (m,1); s = spaugment (A); [L, U, P, Q] = lu (s); x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); x1 = x1(end - n + 1 : end); To find the solution of an overdetermined problem needs an estimate of the residual error R and so it is more complex to formulate a minimum norm solution using the 'spaugment' function. In general the left division operator is more stable and faster than using the 'spaugment' function. See also: *note mldivide: XREFmldivide. Finally, the function 'eigs' can be used to calculate a limited number of eigenvalues and eigenvectors based on a selection criteria and likewise for 'svds' which calculates a limited number of singular values and vectors. -- Function File: D = eigs (A) -- Function File: D = eigs (A, K) -- Function File: D = eigs (A, K, SIGMA) -- Function File: D = eigs (A, K, SIGMA, OPTS) -- Function File: D = eigs (A, B) -- Function File: D = eigs (A, B, K) -- Function File: D = eigs (A, B, K, SIGMA) -- Function File: D = eigs (A, B, K, SIGMA, OPTS) -- Function File: D = eigs (AF, N) -- Function File: D = eigs (AF, N, B) -- Function File: D = eigs (AF, N, K) -- Function File: D = eigs (AF, N, B, K) -- Function File: D = eigs (AF, N, K, SIGMA) -- Function File: D = eigs (AF, N, B, K, SIGMA) -- Function File: D = eigs (AF, N, K, SIGMA, OPTS) -- Function File: D = eigs (AF, N, B, K, SIGMA, OPTS) -- Function File: [V, D] = eigs (A, ...) -- Function File: [V, D] = eigs (AF, N, ...) -- Function File: [V, D, FLAG] = eigs (A, ...) -- Function File: [V, D, FLAG] = eigs (AF, N, ...) Calculate a limited number of eigenvalues and eigenvectors of A, based on a selection criteria. The number of eigenvalues and eigenvectors to calculate is given by K and defaults to 6. By default, 'eigs' solve the equation 'A * v = lambda * v', where 'lambda' is a scalar representing one of the eigenvalues, and 'v' is the corresponding eigenvector. If given the positive definite matrix B then 'eigs' solves the general eigenvalue equation 'A * v = lambda * B * v'. The argument SIGMA determines which eigenvalues are returned. SIGMA can be either a scalar or a string. When SIGMA is a scalar, the K eigenvalues closest to SIGMA are returned. If SIGMA is a string, it must have one of the following values. "lm" Largest Magnitude (default). "sm" Smallest Magnitude. "la" Largest Algebraic (valid only for real symmetric problems). "sa" Smallest Algebraic (valid only for real symmetric problems). "be" Both Ends, with one more from the high-end if K is odd (valid only for real symmetric problems). "lr" Largest Real part (valid only for complex or unsymmetric problems). "sr" Smallest Real part (valid only for complex or unsymmetric problems). "li" Largest Imaginary part (valid only for complex or unsymmetric problems). "si" Smallest Imaginary part (valid only for complex or unsymmetric problems). If OPTS is given, it is a structure defining possible options that 'eigs' should use. The fields of the OPTS structure are: 'issym' If AF is given, then flags whether the function AF defines a symmetric problem. It is ignored if A is given. The default is false. 'isreal' If AF is given, then flags whether the function AF defines a real problem. It is ignored if A is given. The default is true. 'tol' Defines the required convergence tolerance, calculated as 'tol * norm (A)'. The default is 'eps'. 'maxit' The maximum number of iterations. The default is 300. 'p' The number of Lanzcos basis vectors to use. More vectors will result in faster convergence, but a greater use of memory. The optimal value of 'p' is problem dependent and should be in the range K to N. The default value is '2 * K'. 'v0' The starting vector for the algorithm. An initial vector close to the final vector will speed up convergence. The default is for ARPACK to randomly generate a starting vector. If specified, 'v0' must be an N-by-1 vector where 'N = rows (A)' 'disp' The level of diagnostic printout (0|1|2). If 'disp' is 0 then diagnostics are disabled. The default value is 0. 'cholB' Flag if 'chol (B)' is passed rather than B. The default is false. 'permB' The permutation vector of the Cholesky factorization of B if 'cholB' is true. That is 'chol (B(permB, permB))'. The default is '1:N'. It is also possible to represent A by a function denoted AF. AF must be followed by a scalar argument N defining the length of the vector argument accepted by AF. AF can be a function handle, an inline function, or a string. When AF is a string it holds the name of the function to use. AF is a function of the form 'y = af (x)' where the required return value of AF is determined by the value of SIGMA. The four possible forms are 'A * x' if SIGMA is not given or is a string other than "sm". 'A \ x' if SIGMA is 0 or "sm". '(A - sigma * I) \ x' for the standard eigenvalue problem, where 'I' is the identity matrix of the same size as A. '(A - sigma * B) \ x' for the general eigenvalue problem. The return arguments of 'eigs' depend on the number of return arguments requested. With a single return argument, a vector D of length K is returned containing the K eigenvalues that have been found. With two return arguments, V is a N-by-K matrix whose columns are the K eigenvectors corresponding to the returned eigenvalues. The eigenvalues themselves are returned in D in the form of a N-by-K matrix, where the elements on the diagonal are the eigenvalues. Given a third return argument FLAG, 'eigs' returns the status of the convergence. If FLAG is 0 then all eigenvalues have converged. Any other value indicates a failure to converge. This function is based on the ARPACK package, written by R. Lehoucq, K. Maschhoff, D. Sorensen, and C. Yang. For more information see . See also: *note eig: XREFeig, *note svds: XREFsvds. -- Function File: S = svds (A) -- Function File: S = svds (A, K) -- Function File: S = svds (A, K, SIGMA) -- Function File: S = svds (A, K, SIGMA, OPTS) -- Function File: [U, S, V] = svds (...) -- Function File: [U, S, V, FLAG] = svds (...) Find a few singular values of the matrix A. The singular values are calculated using [M, N] = size (A); S = eigs ([sparse(M, M), A; A', sparse(N, N)]) The eigenvalues returned by 'eigs' correspond to the singular values of A. The number of singular values to calculate is given by K and defaults to 6. The argument SIGMA specifies which singular values to find. When SIGMA is the string 'L', the default, the largest singular values of A are found. Otherwise, SIGMA must be a real scalar and the singular values closest to SIGMA are found. As a corollary, 'SIGMA = 0' finds the smallest singular values. Note that for relatively small values of SIGMA, there is a chance that the requested number of singular values will not be found. In that case SIGMA should be increased. OPTS is a structure defining options that 'svds' will pass to 'eigs'. The possible fields of this structure are documented in 'eigs'. By default, 'svds' sets the following three fields: 'tol' The required convergence tolerance for the singular values. The default value is 1e-10. 'eigs' is passed 'TOL / sqrt(2)'. 'maxit' The maximum number of iterations. The default is 300. 'disp' The level of diagnostic printout (0|1|2). If 'disp' is 0 then diagnostics are disabled. The default value is 0. If more than one output is requested then 'svds' will return an approximation of the singular value decomposition of A A_approx = U*S*V' where A_approx is a matrix of size A but only rank K. FLAG returns 0 if the algorithm has succesfully converged, and 1 otherwise. The test for convergence is norm (A*V - U*S, 1) <= TOL * norm (A, 1) 'svds' is best for finding only a few singular values from a large sparse matrix. Otherwise, 'svd (full (A))' will likely be more efficient. See also: *note svd: XREFsvd, *note eigs: XREFeigs. ---------- Footnotes ---------- (1) The CHOLMOD, UMFPACK and CXSPARSE packages were written by Tim Davis and are available at  File: octave.info, Node: Iterative Techniques, Next: Real Life Example, Prev: Sparse Linear Algebra, Up: Sparse Matrices 22.3 Iterative Techniques Applied to Sparse Matrices ==================================================== The left division '\' and right division '/' operators, discussed in the previous section, use direct solvers to resolve a linear equation of the form 'X = A \ B' or 'X = B / A'. Octave also includes a number of functions to solve sparse linear equations using iterative techniques. -- Function File: X = pcg (A, B, TOL, MAXIT, M1, M2, X0, ...) -- Function File: [X, FLAG, RELRES, ITER, RESVEC, EIGEST] = pcg (...) Solve the linear system of equations 'A * X = B' by means of the Preconditioned Conjugate Gradient iterative method. The input arguments are * A can be either a square (preferably sparse) matrix or a function handle, inline function or string containing the name of a function which computes 'A * X'. In principle, A should be symmetric and positive definite; if 'pcg' finds A not to be positive definite, a warning is printed and the FLAG output will be set. * B is the right-hand side vector. * TOL is the required relative tolerance for the residual error, 'B - A * X'. The iteration stops if 'norm (B - A * X)' <= TOL * norm (B). If TOL is omitted or empty then a tolerance of 1e-6 is used. * MAXIT is the maximum allowable number of iterations; if MAXIT is omitted or empty then a value of 20 is used. * M = M1 * M2 is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by 'pcg' 'P * X = M \ B', with 'P = M \ A'. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method. Instead of matrices M1 and M2, the user may pass two functions which return the results of applying the inverse of M1 and M2 to a vector (usually this is the preferred way of using the preconditioner). If M1 is omitted or empty '[]' then no preconditioning is applied. If M2 is omitted, M = M1 will be used as a preconditioner. * X0 is the initial guess. If X0 is omitted or empty then the function sets X0 to a zero vector by default. The arguments which follow X0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to 'pcg'. See the examples below for further details. The output arguments are * X is the computed approximation to the solution of 'A * X = B'. * FLAG reports on the convergence. A value of 0 means the solution converged and the tolerance criterion given by TOL is satisfied. A value of 1 means that the MAXIT limit for the iteration count was reached. A value of 3 indicates that the (preconditioned) matrix was found not to be positive definite. * RELRES is the ratio of the final residual to its initial value, measured in the Euclidean norm. * ITER is the actual number of iterations performed. * RESVEC describes the convergence history of the method. 'RESVEC(i,1)' is the Euclidean norm of the residual, and 'RESVEC(i,2)' is the preconditioned residual norm, after the (I-1)-th iteration, 'I = 1, 2, ..., ITER+1'. The preconditioned residual norm is defined as 'norm (R) ^ 2 = R' * (M \ R)' where 'R = B - A * X', see also the description of M. If EIGEST is not required, only 'RESVEC(:,1)' is returned. * EIGEST returns the estimate for the smallest 'EIGEST(1)' and largest 'EIGEST(2)' eigenvalues of the preconditioned matrix 'P = M \ A'. In particular, if no preconditioning is used, the estimates for the extreme eigenvalues of A are returned. 'EIGEST(1)' is an overestimate and 'EIGEST(2)' is an underestimate, so that 'EIGEST(2) / EIGEST(1)' is a lower bound for 'cond (P, 2)', which nevertheless in the limit should theoretically be equal to the actual value of the condition number. The method which computes EIGEST works only for symmetric positive definite A and M, and the user is responsible for verifying this assumption. Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A) n = 10; A = diag (sparse (1:n)); b = rand (n, 1); [l, u, p] = ilu (A, struct ("droptol", 1.e-3)); EXAMPLE 1: Simplest use of 'pcg' x = pcg (A, b) EXAMPLE 2: 'pcg' with a function which computes 'A * X' function y = apply_a (x) y = [1:N]' .* x; endfunction x = pcg ("apply_a", b) EXAMPLE 3: 'pcg' with a preconditioner: L * U x = pcg (A, b, 1.e-6, 500, l*u) EXAMPLE 4: 'pcg' with a preconditioner: L * U. Faster than EXAMPLE 3 since lower and upper triangular matrices are easier to invert x = pcg (A, b, 1.e-6, 500, l, u) EXAMPLE 5: Preconditioned iteration, with full diagnostics. The preconditioner (quite strange, because even the original matrix A is trivial) is defined as a function function y = apply_m (x) k = floor (length (x) - 2); y = x; y(1:k) = x(1:k) ./ [1:k]'; endfunction [x, flag, relres, iter, resvec, eigest] = ... pcg (A, b, [], [], "apply_m"); semilogy (1:iter+1, resvec); EXAMPLE 6: Finally, a preconditioner which depends on a parameter K. function y = apply_M (x, varargin) K = varargin{1}; y = x; y(1:K) = x(1:K) ./ [1:K]'; endfunction [x, flag, relres, iter, resvec, eigest] = ... pcg (A, b, [], [], "apply_m", [], [], 3) References: 1. C.T. Kelley, 'Iterative Methods for Linear and Nonlinear Equations', SIAM, 1995. (the base PCG algorithm) 2. Y. Saad, 'Iterative Methods for Sparse Linear Systems', PWS 1996. (condition number estimate from PCG) Revised version of this book is available online at See also: *note sparse: XREFsparse, *note pcr: XREFpcr. -- Function File: X = pcr (A, B, TOL, MAXIT, M, X0, ...) -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = pcr (...) Solve the linear system of equations 'A * X = B' by means of the Preconditioned Conjugate Residuals iterative method. The input arguments are * A can be either a square (preferably sparse) matrix or a function handle, inline function or string containing the name of a function which computes 'A * X'. In principle A should be symmetric and non-singular; if 'pcr' finds A to be numerically singular, you will get a warning message and the FLAG output parameter will be set. * B is the right hand side vector. * TOL is the required relative tolerance for the residual error, 'B - A * X'. The iteration stops if 'norm (B - A * X) <= TOL * norm (B - A * X0)'. If TOL is empty or is omitted, the function sets 'TOL = 1e-6' by default. * MAXIT is the maximum allowable number of iterations; if '[]' is supplied for 'maxit', or 'pcr' has less arguments, a default value equal to 20 is used. * M is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by 'pcr' 'P * X = M \ B', with 'P = M \ A'. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method. Instead of matrix M, the user may pass a function which returns the results of applying the inverse of M to a vector (usually this is the preferred way of using the preconditioner). If '[]' is supplied for M, or M is omitted, no preconditioning is applied. * X0 is the initial guess. If X0 is empty or omitted, the function sets X0 to a zero vector by default. The arguments which follow X0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to 'pcr'. See the examples below for further details. The output arguments are * X is the computed approximation to the solution of 'A * X = B'. * FLAG reports on the convergence. 'FLAG = 0' means the solution converged and the tolerance criterion given by TOL is satisfied. 'FLAG = 1' means that the MAXIT limit for the iteration count was reached. 'FLAG = 3' reports a 'pcr' breakdown, see [1] for details. * RELRES is the ratio of the final residual to its initial value, measured in the Euclidean norm. * ITER is the actual number of iterations performed. * RESVEC describes the convergence history of the method, so that 'RESVEC (i)' contains the Euclidean norms of the residual after the (I-1)-th iteration, 'I = 1,2, ..., ITER+1'. Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A) n = 10; A = sparse (diag (1:n)); b = rand (N, 1); EXAMPLE 1: Simplest use of 'pcr' x = pcr (A, b) EXAMPLE 2: 'pcr' with a function which computes 'A * X'. function y = apply_a (x) y = [1:10]' .* x; endfunction x = pcr ("apply_a", b) EXAMPLE 3: Preconditioned iteration, with full diagnostics. The preconditioner (quite strange, because even the original matrix A is trivial) is defined as a function function y = apply_m (x) k = floor (length (x) - 2); y = x; y(1:k) = x(1:k) ./ [1:k]'; endfunction [x, flag, relres, iter, resvec] = ... pcr (A, b, [], [], "apply_m") semilogy ([1:iter+1], resvec); EXAMPLE 4: Finally, a preconditioner which depends on a parameter K. function y = apply_m (x, varargin) k = varargin{1}; y = x; y(1:k) = x(1:k) ./ [1:k]'; endfunction [x, flag, relres, iter, resvec] = ... pcr (A, b, [], [], "apply_m"', [], 3) References: [1] W. Hackbusch, 'Iterative Solution of Large Sparse Systems of Equations', section 9.5.4; Springer, 1994 See also: *note sparse: XREFsparse, *note pcg: XREFpcg. The speed with which an iterative solver converges to a solution can be accelerated with the use of a pre-conditioning matrix M. In this case the linear equation 'M^-1 * X = M^-1 * A \ B' is solved instead. Typical pre-conditioning matrices are partial factorizations of the original matrix. -- Function File: L = ichol (A) -- Function File: L = ichol (A, OPTS) Compute the incomplete Cholesky factorization of the sparse square matrix A. By default, 'ichol' uses only the lower triangle of A and produces a lower triangular factor L such that L*L' approximates A. The factor given by this routine may be useful as a preconditioner for a system of linear equations being solved by iterative methods such as PCG (Preconditioned Conjugate Gradient). The factorization may be modified by passing options in a structure OPTS. The option name is a field of the structure and the setting is the value of field. Names and specifiers are case sensitive. type Type of factorization. "nofill" (default) Incomplete Cholesky factorization with no fill-in (IC(0)). "ict" Incomplete Cholesky factorization with threshold dropping (ICT). diagcomp A non-negative scalar ALPHA for incomplete Cholesky factorization of 'A + ALPHA * diag (diag (A))' instead of A. This can be useful when A is not positive definite. The default value is 0. droptol A non-negative scalar specifying the drop tolerance for factorization if performing ICT. The default value is 0 which produces the complete Cholesky factorization. Non-diagonal entries of L are set to 0 unless 'abs (L(i,j)) >= droptol * norm (A(j:end, j), 1)'. michol Modified incomplete Cholesky factorization: "off" (default) Row and column sums are not necessarily preserved. "on" The diagonal of L is modified so that row (and column) sums are preserved even when elements have been dropped during the factorization. The relationship preserved is: 'A * e = L * L' * e', where e is a vector of ones. shape "lower" (default) Use only the lower triangle of A and return a lower triangular factor L such that L*L' approximates A. "upper" Use only the upper triangle of A and return an upper triangular factor U such that 'U'*U' approximates A. EXAMPLES The following problem demonstrates how to factorize a sample symmetric positive definite matrix with the full Cholesky decomposition and with the incomplete one. A = [ 0.37, -0.05, -0.05, -0.07; -0.05, 0.116, 0.0, -0.05; -0.05, 0.0, 0.116, -0.05; -0.07, -0.05, -0.05, 0.202]; A = sparse (A); nnz (tril (A)) ans = 9 L = chol (A, "lower"); nnz (L) ans = 10 norm (A - L * L', "fro") / norm (A, "fro") ans = 1.1993e-16 opts.type = "nofill"; L = ichol (A, opts); nnz (L) ans = 9 norm (A - L * L', "fro") / norm (A, "fro") ans = 0.019736 Another example for decomposition is a finite difference matrix used to solve a boundary value problem on the unit square. nx = 400; ny = 200; hx = 1 / (nx + 1); hy = 1 / (ny + 1); Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2); Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2); A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy); nnz (tril (A)) ans = 239400 opts.type = "nofill"; L = ichol (A, opts); nnz (tril (A)) ans = 239400 norm (A - L * L', "fro") / norm (A, "fro") ans = 0.062327 References for implemented algorithms: [1] Y. Saad. "Preconditioning Techniques." 'Iterative Methods for Sparse Linear Systems', PWS Publishing Company, 1996. [2] M. Jones, P. Plassmann: 'An Improved Incomplete Cholesky Factorization', 1992. See also: *note chol: XREFchol, *note ilu: XREFilu, *note pcg: XREFpcg. -- Function File: ilu (A) -- Function File: ilu (A, OPTS) -- Function File: [L, U] = ilu (...) -- Function File: [L, U, P] = ilu (...) Compute the incomplete LU factorization of the sparse square matrix A. 'ilu' returns a unit lower triangular matrix L, an upper triangular matrix U, and optionally a permutation matrix P, such that 'L*U' approximates 'P*A'. The factors given by this routine may be useful as preconditioners for a system of linear equations being solved by iterative methods such as BICG (BiConjugate Gradients) or GMRES (Generalized Minimum Residual Method). The factorization may be modified by passing options in a structure OPTS. The option name is a field of the structure and the setting is the value of field. Names and specifiers are case sensitive. 'type' Type of factorization. "nofill" ILU factorization with no fill-in (ILU(0)). Additional supported options: 'milu'. "crout" Crout version of ILU factorization (ILUC). Additional supported options: 'milu', 'droptol'. "ilutp" (default) ILU factorization with threshold and pivoting. Additional supported options: 'milu', 'droptol', 'udiag', 'thresh'. 'droptol' A non-negative scalar specifying the drop tolerance for factorization. The default value is 0 which produces the complete LU factorization. Non-diagonal entries of U are set to 0 unless 'abs (U(i,j)) >= droptol * norm (A(:,j))'. Non-diagonal entries of L are set to 0 unless 'abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j)'. 'milu' Modified incomplete LU factorization: "row" Row-sum modified incomplete LU factorization. The factorization preserves row sums: 'A * e = L * U * e', where e is a vector of ones. "col" Column-sum modified incomplete LU factorization. The factorization preserves column sums: 'e' * A = e' * L * U'. "off" (default) Row and column sums are not necessarily preserved. 'udiag' If true, any zeros on the diagonal of the upper triangular factor are replaced by the local drop tolerance 'droptol * norm (A(:,j))/U(j,j)'. The default is false. 'thresh' Pivot threshold for factorization. It can range between 0 (diagonal pivoting) and 1 (default), where the maximum magnitude entry in the column is chosen to be the pivot. If 'ilu' is called with just one output, the returned matrix is 'L + U - speye (size (A))', where L is unit lower triangular and U is upper triangular. With two outputs, 'ilu' returns a unit lower triangular matrix L and an upper triangular matrix U. For OPTS.type == "ilutp", one of the factors is permuted based on the value of OPTS.milu. When OPTS.milu == "row", U is a column permuted upper triangular factor. Otherwise, L is a row-permuted unit lower triangular factor. If there are three named outputs and OPTS.milu != "row", P is returned such that L and U are incomplete factors of 'P*A'. When OPTS.milu == "row", P is returned such that L and U are incomplete factors of 'A*P'. EXAMPLES A = gallery ("neumann", 1600) + speye (1600); opts.type = "nofill"; nnz (A) ans = 7840 nnz (lu (A)) ans = 126478 nnz (ilu (A, opts)) ans = 7840 This shows that A has 7,840 nonzeros, the complete LU factorization has 126,478 nonzeros, and the incomplete LU factorization, with 0 level of fill-in, has 7,840 nonzeros, the same amount as A. Taken from: http://www.mathworks.com/help/matlab/ref/ilu.html A = gallery ("wathen", 10, 10); b = sum (A, 2); tol = 1e-8; maxit = 50; opts.type = "crout"; opts.droptol = 1e-4; [L, U] = ilu (A, opts); x = bicg (A, b, tol, maxit, L, U); norm (A * x - b, inf) This example uses ILU as preconditioner for a random FEM-Matrix, which has a large condition number. Without L and U BICG would not converge. See also: *note lu: XREFlu, *note ichol: XREFichol, *note bicg: XREFbicg, *note gmres: XREFgmres.  File: octave.info, Node: Real Life Example, Prev: Iterative Techniques, Up: Sparse Matrices 22.4 Real Life Example using Sparse Matrices ============================================ A common application for sparse matrices is in the solution of Finite Element Models. Finite element models allow numerical solution of partial differential equations that do not have closed form solutions, typically because of the complex shape of the domain. In order to motivate this application, we consider the boundary value Laplace equation. This system can model scalar potential fields, such as heat or electrical potential. Given a medium Omega with boundary dOmega. At all points on the dOmega the boundary conditions are known, and we wish to calculate the potential in Omega. Boundary conditions may specify the potential (Dirichlet boundary condition), its normal derivative across the boundary (Neumann boundary condition), or a weighted sum of the potential and its derivative (Cauchy boundary condition). In a thermal model, we want to calculate the temperature in Omega and know the boundary temperature (Dirichlet condition) or heat flux (from which we can calculate the Neumann condition by dividing by the thermal conductivity at the boundary). Similarly, in an electrical model, we want to calculate the voltage in Omega and know the boundary voltage (Dirichlet) or current (Neumann condition after diving by the electrical conductivity). In an electrical model, it is common for much of the boundary to be electrically isolated; this is a Neumann boundary condition with the current equal to zero. The simplest finite element models will divide Omega into simplexes (triangles in 2D, pyramids in 3D). The following example creates a simple rectangular 2-D electrically conductive medium with 10 V and 20 V imposed on opposite sides (Dirichlet boundary conditions). All other edges are electrically isolated. node_y = [1;1.2;1.5;1.8;2]*ones(1,11); node_x = ones(5,1)*[1,1.05,1.1,1.2, ... 1.3,1.5,1.7,1.8,1.9,1.95,2]; nodes = [node_x(:), node_y(:)]; [h,w] = size (node_x); elems = []; for idx = 1:w-1 widx = (idx-1)*h; elems = [elems; ... widx+[(1:h-1);(2:h);h+(1:h-1)]'; ... widx+[(2:h);h+(2:h);h+(1:h-1)]' ]; endfor E = size (elems,1); # No. of simplices N = size (nodes,1); # No. of vertices D = size (elems,2); # dimensions+1 This creates a N-by-2 matrix 'nodes' and a E-by-3 matrix 'elems' with values, which define finite element triangles: nodes(1:7,:)' 1.00 1.00 1.00 1.00 1.00 1.05 1.05 ... 1.00 1.20 1.50 1.80 2.00 1.00 1.20 ... elems(1:7,:)' 1 2 3 4 2 3 4 ... 2 3 4 5 7 8 9 ... 6 7 8 9 6 7 8 ... Using a first order FEM, we approximate the electrical conductivity distribution in Omega as constant on each simplex (represented by the vector 'conductivity'). Based on the finite element geometry, we first calculate a system (or stiffness) matrix for each simplex (represented as 3-by-3 elements on the diagonal of the element-wise system matrix 'SE'). Based on 'SE' and a N-by-DE connectivity matrix 'C', representing the connections between simplices and vertices, the global connectivity matrix 'S' is calculated. ## Element conductivity conductivity = [1*ones(1,16), ... 2*ones(1,48), 1*ones(1,16)]; ## Connectivity matrix C = sparse ((1:D*E), reshape (elems', ... D*E, 1), 1, D*E, N); ## Calculate system matrix Siidx = floor ([0:D*E-1]'/D) * D * ... ones(1,D) + ones(D*E,1)*(1:D) ; Sjidx = [1:D*E]'*ones (1,D); Sdata = zeros (D*E,D); dfact = factorial (D-1); for j = 1:E a = inv ([ones(D,1), ... nodes(elems(j,:), :)]); const = conductivity(j) * 2 / ... dfact / abs (det (a)); Sdata(D*(j-1)+(1:D),:) = const * ... a(2:D,:)' * a(2:D,:); endfor ## Element-wise system matrix SE = sparse(Siidx,Sjidx,Sdata); ## Global system matrix S = C'* SE *C; The system matrix acts like the conductivity 'S' in Ohm's law 'S * V = I'. Based on the Dirichlet and Neumann boundary conditions, we are able to solve for the voltages at each vertex 'V'. ## Dirichlet boundary conditions D_nodes = [1:5, 51:55]; D_value = [10*ones(1,5), 20*ones(1,5)]; V = zeros (N,1); V(D_nodes) = D_value; idx = 1:N; # vertices without Dirichlet # boundary condns idx(D_nodes) = []; ## Neumann boundary conditions. Note that ## N_value must be normalized by the ## boundary length and element conductivity N_nodes = []; N_value = []; Q = zeros (N,1); Q(N_nodes) = N_value; V(idx) = S(idx,idx) \ ( Q(idx) - ... S(idx,D_nodes) * V(D_nodes)); Finally, in order to display the solution, we show each solved voltage value in the z-axis for each simplex vertex. elemx = elems(:,[1,2,3,1])'; xelems = reshape (nodes(elemx, 1), 4, E); yelems = reshape (nodes(elemx, 2), 4, E); velems = reshape (V(elemx), 4, E); plot3 (xelems,yelems,velems,"k"); print "grid.eps";  File: octave.info, Node: Numerical Integration, Next: Differential Equations, Prev: Sparse Matrices, Up: Top 23 Numerical Integration ************************ Octave comes with several built-in functions for computing the integral of a function numerically (termed quadrature). These functions all solve 1-dimensional integration problems. * Menu: * Functions of One Variable:: * Orthogonal Collocation:: * Functions of Multiple Variables::  File: octave.info, Node: Functions of One Variable, Next: Orthogonal Collocation, Up: Numerical Integration 23.1 Functions of One Variable ============================== Octave supports five different algorithms for computing the integral of a function f over the interval from a to b. These are 'quad' Numerical integration based on Gaussian quadrature. 'quadv' Numerical integration using an adaptive vectorized Simpson's rule. 'quadl' Numerical integration using an adaptive Lobatto rule. 'quadgk' Numerical integration using an adaptive Gauss-Konrod rule. 'quadcc' Numerical integration using adaptive Clenshaw-Curtis rules. 'trapz, cumtrapz' Numerical integration of data using the trapezoidal method. The best quadrature algorithm to use depends on the integrand. If you have empirical data, rather than a function, the choice is 'trapz' or 'cumtrapz'. If you are uncertain about the characteristics of the integrand, 'quadcc' will be the most robust as it can handle discontinuities, singularities, oscillatory functions, and infinite intervals. When the integrand is smooth 'quadgk' may be the fastest of the algorithms. Function Characteristics ---------------------------------------------------------------------------- quad Low accuracy with nonsmooth integrands quadv Medium accuracy with smooth integrands quadl Medium accuracy with smooth integrands. Slower than quadgk. quadgk Medium accuracy (1e^{-6}-1e^{-9}) with smooth integrands. Handles oscillatory functions and infinite bounds quadcc Low to High accuracy with nonsmooth/smooth integrands Handles oscillatory functions, singularities, and infinite bounds Here is an example of using 'quad' to integrate the function F(X) = X * sin (1/X) * sqrt (abs (1 - X)) from X = 0 to X = 3. This is a fairly difficult integration (plot the function over the range of integration to see why). The first step is to define the function: function y = f (x) y = x .* sin (1./x) .* sqrt (abs (1 - x)); endfunction Note the use of the 'dot' forms of the operators. This is not necessary for the 'quad' integrator, but is required by the other integrators. In any case, it makes it much easier to generate a set of points for plotting because it is possible to call the function with a vector argument to produce a vector result. The second step is to call quad with the limits of integration: [q, ier, nfun, err] = quad ("f", 0, 3) => 1.9819 => 1 => 5061 => 1.1522e-07 Although 'quad' returns a nonzero value for IER, the result is reasonably accurate (to see why, examine what happens to the result if you move the lower bound to 0.1, then 0.01, then 0.001, etc.). The function "f" can be the string name of a function, a function handle, or an inline function. These options make it quite easy to do integration without having to fully define a function in an m-file. For example: # Verify integral (x^3) = x^4/4 f = inline ("x.^3"); quadgk (f, 0, 1) => 0.25000 # Verify gamma function = (n-1)! for n = 4 f = @(x) x.^3 .* exp (-x); quadcc (f, 0, Inf) => 6.0000 -- Built-in Function: Q = quad (F, A, B) -- Built-in Function: Q = quad (F, A, B, TOL) -- Built-in Function: Q = quad (F, A, B, TOL, SING) -- Built-in Function: [Q, IER, NFUN, ERR] = quad (...) Numerically evaluate the integral of F from A to B using Fortran routines from QUADPACK. F is a function handle, inline function, or a string containing the name of the function to evaluate. The function must have the form 'y = f (x)' where Y and X are scalars. A and B are the lower and upper limits of integration. Either or both may be infinite. The optional argument TOL is a vector that specifies the desired accuracy of the result. The first element of the vector is the desired absolute tolerance, and the second element is the desired relative tolerance. To choose a relative test only, set the absolute tolerance to zero. To choose an absolute test only, set the relative tolerance to zero. Both tolerances default to 'sqrt (eps)' or approximately 1.5e^{-8}. The optional argument SING is a vector of values at which the integrand is known to be singular. The result of the integration is returned in Q. IER contains an integer error code (0 indicates a successful integration). NFUN indicates the number of function evaluations that were made. ERR contains an estimate of the error in the solution. The function 'quad_options' can set other optional parameters for 'quad'. Note: because 'quad' is written in Fortran it cannot be called recursively. This prevents its use in integrating over more than one variable by routines 'dblquad' and 'triplequad'. See also: *note quad_options: XREFquad_options, *note quadv: XREFquadv, *note quadl: XREFquadl, *note quadgk: XREFquadgk, *note quadcc: XREFquadcc, *note trapz: XREFtrapz, *note dblquad: XREFdblquad, *note triplequad: XREFtriplequad. -- Built-in Function: quad_options () -- Built-in Function: val = quad_options (OPT) -- Built-in Function: quad_options (OPT, VAL) Query or set options for the function 'quad'. When called with no arguments, the names of all available options and their current values are displayed. Given one argument, return the value of the option OPT. When called with two arguments, 'quad_options' sets the option OPT to value VAL. Options include '"absolute tolerance"' Absolute tolerance; may be zero for pure relative error test. '"relative tolerance"' Non-negative relative tolerance. If the absolute tolerance is zero, the relative tolerance must be greater than or equal to 'max (50*eps, 0.5e-28)'. '"single precision absolute tolerance"' Absolute tolerance for single precision; may be zero for pure relative error test. '"single precision relative tolerance"' Non-negative relative tolerance for single precision. If the absolute tolerance is zero, the relative tolerance must be greater than or equal to 'max (50*eps, 0.5e-28)'. -- Function File: Q = quadv (F, A, B) -- Function File: Q = quadv (F, A, B, TOL) -- Function File: Q = quadv (F, A, B, TOL, TRACE) -- Function File: Q = quadv (F, A, B, TOL, TRACE, P1, P2, ...) -- Function File: [Q, NFUN] = quadv (...) Numerically evaluate the integral of F from A to B using an adaptive Simpson's rule. F is a function handle, inline function, or string containing the name of the function to evaluate. 'quadv' is a vectorized version of 'quad' and the function defined by F must accept a scalar or vector as input and return a scalar, vector, or array as output. A and B are the lower and upper limits of integration. Both limits must be finite. The optional argument TOL defines the tolerance used to stop the adaptation procedure. The default value is 1e-6. The algorithm used by 'quadv' involves recursively subdividing the integration interval and applying Simpson's rule on each subinterval. If TRACE is true then after computing each of these partial integrals display: (1) the total number of function evaluations, (2) the left end of the subinterval, (3) the length of the subinterval, (4) the approximation of the integral over the subinterval. Additional arguments P1, etc., are passed directly to the function F. To use default values for TOL and TRACE, one may pass empty matrices ([]). The result of the integration is returned in Q NFUN indicates the number of function evaluations that were made. Note: 'quadv' is written in Octave's scripting language and can be used recursively in 'dblquad' and 'triplequad', unlike the 'quad' function. See also: *note quad: XREFquad, *note quadl: XREFquadl, *note quadgk: XREFquadgk, *note quadcc: XREFquadcc, *note trapz: XREFtrapz, *note dblquad: XREFdblquad, *note triplequad: XREFtriplequad. -- Function File: Q = quadl (F, A, B) -- Function File: Q = quadl (F, A, B, TOL) -- Function File: Q = quadl (F, A, B, TOL, TRACE) -- Function File: Q = quadl (F, A, B, TOL, TRACE, P1, P2, ...) Numerically evaluate the integral of F from A to B using an adaptive Lobatto rule. F is a function handle, inline function, or string containing the name of the function to evaluate. The function F must be vectorized and return a vector of output values when given a vector of input values. A and B are the lower and upper limits of integration. Both limits must be finite. The optional argument TOL defines the relative tolerance with which to perform the integration. The default value is 'eps'. The algorithm used by 'quadl' involves recursively subdividing the integration interval. If TRACE is defined then for each subinterval display: (1) the left end of the subinterval, (2) the length of the subinterval, (3) the approximation of the integral over the subinterval. Additional arguments P1, etc., are passed directly to the function F. To use default values for TOL and TRACE, one may pass empty matrices ([]). Reference: W. Gander and W. Gautschi, 'Adaptive Quadrature - Revisited', BIT Vol. 40, No. 1, March 2000, pp. 84-101. See also: *note quad: XREFquad, *note quadv: XREFquadv, *note quadgk: XREFquadgk, *note quadcc: XREFquadcc, *note trapz: XREFtrapz, *note dblquad: XREFdblquad, *note triplequad: XREFtriplequad. -- Function File: Q = quadgk (F, A, B) -- Function File: Q = quadgk (F, A, B, ABSTOL) -- Function File: Q = quadgk (F, A, B, ABSTOL, TRACE) -- Function File: Q = quadgk (F, A, B, PROP, VAL, ...) -- Function File: [Q, ERR] = quadgk (...) Numerically evaluate the integral of F from A to B using adaptive Gauss-Konrod quadrature. F is a function handle, inline function, or string containing the name of the function to evaluate. The function F must be vectorized and return a vector of output values when given a vector of input values. A and B are the lower and upper limits of integration. Either or both limits may be infinite or contain weak end singularities. Variable transformation will be used to treat any infinite intervals and weaken the singularities. For example: quadgk (@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf) Note that the formulation of the integrand uses the element-by-element operator './' and all user functions to 'quadgk' should do the same. The optional argument TOL defines the absolute tolerance used to stop the integration procedure. The default value is 1e-10. The algorithm used by 'quadgk' involves subdividing the integration interval and evaluating each subinterval. If TRACE is true then after computing each of these partial integrals display: (1) the number of subintervals at this step, (2) the current estimate of the error ERR, (3) the current estimate for the integral Q. Alternatively, properties of 'quadgk' can be passed to the function as pairs "PROP", VAL. Valid properties are 'AbsTol' Define the absolute error tolerance for the quadrature. The default absolute tolerance is 1e-10. 'RelTol' Define the relative error tolerance for the quadrature. The default relative tolerance is 1e-5. 'MaxIntervalCount' 'quadgk' initially subdivides the interval on which to perform the quadrature into 10 intervals. Subintervals that have an unacceptable error are subdivided and re-evaluated. If the number of subintervals exceeds 650 subintervals at any point then a poor convergence is signaled and the current estimate of the integral is returned. The property "MaxIntervalCount" can be used to alter the number of subintervals that can exist before exiting. 'WayPoints' Discontinuities in the first derivative of the function to integrate can be flagged with the "WayPoints" property. This forces the ends of a subinterval to fall on the breakpoints of the function and can result in significantly improved estimation of the error in the integral, faster computation, or both. For example, quadgk (@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1) signals the breakpoint in the integrand at 'X = 1'. 'Trace' If logically true 'quadgk' prints information on the convergence of the quadrature at each iteration. If any of A, B, or WAYPOINTS is complex then the quadrature is treated as a contour integral along a piecewise continuous path defined by the above. In this case the integral is assumed to have no edge singularities. For example, quadgk (@(z) log (z), 1+1i, 1+1i, "WayPoints", [1-1i, -1,-1i, -1+1i]) integrates 'log (z)' along the square defined by '[1+1i, 1-1i, -1-1i, -1+1i]'. The result of the integration is returned in Q. ERR is an approximate bound on the error in the integral 'abs (Q - I)', where I is the exact value of the integral. Reference: L.F. Shampine, '"Vectorized adaptive quadrature in MATLAB"', Journal of Computational and Applied Mathematics, pp. 131-140, Vol 211, Issue 2, Feb 2008. See also: *note quad: XREFquad, *note quadv: XREFquadv, *note quadl: XREFquadl, *note quadcc: XREFquadcc, *note trapz: XREFtrapz, *note dblquad: XREFdblquad, *note triplequad: XREFtriplequad. -- Function File: Q = quadcc (F, A, B) -- Function File: Q = quadcc (F, A, B, TOL) -- Function File: Q = quadcc (F, A, B, TOL, SING) -- Function File: [Q, ERR, NR_POINTS] = quadcc (...) Numerically evaluate the integral of F from A to B using doubly-adaptive Clenshaw-Curtis quadrature. F is a function handle, inline function, or string containing the name of the function to evaluate. The function F must be vectorized and must return a vector of output values if given a vector of input values. For example, f = @(x) x .* sin (1./x) .* sqrt (abs (1 - x)); which uses the element-by-element "dot" form for all operators. A and B are the lower and upper limits of integration. Either or both limits may be infinite. 'quadcc' handles an inifinite limit by substituting the variable of integration with 'x = tan (pi/2*u)'. The optional argument TOL defines the relative tolerance used to stop the integration procedure. The default value is 1e^{-6}. The optional argument SING contains a list of points where the integrand has known singularities, or discontinuities in any of its derivatives, inside the integration interval. For the example above, which has a discontinuity at x=1, the call to 'quadcc' would be as follows int = quadcc (f, a, b, 1.0e-6, [ 1 ]); The result of the integration is returned in Q. ERR is an estimate of the absolute integration error. NR_POINTS is the number of points at which the integrand was evaluated. If the adaptive integration did not converge, the value of ERR will be larger than the requested tolerance. Therefore, it is recommended to verify this value for difficult integrands. 'quadcc' is capable of dealing with non-numeric values of the integrand such as 'NaN' or 'Inf'. If the integral diverges, and 'quadcc' detects this, then a warning is issued and 'Inf' or '-Inf' is returned. Note: 'quadcc' is a general purpose quadrature algorithm and, as such, may be less efficient for a smooth or otherwise well-behaved integrand than other methods such as 'quadgk'. The algorithm uses Clenshaw-Curtis quadrature rules of increasing degree in each interval and bisects the interval if either the function does not appear to be smooth or a rule of maximum degree has been reached. The error estimate is computed from the L2-norm of the difference between two successive interpolations of the integrand over the nodes of the respective quadrature rules. Reference: P. Gonnet, 'Increasing the Reliability of Adaptive Quadrature Using Explicit Interpolants', ACM Transactions on Mathematical Software, Vol. 37, Issue 3, Article No. 3, 2010. See also: *note quad: XREFquad, *note quadv: XREFquadv, *note quadl: XREFquadl, *note quadgk: XREFquadgk, *note trapz: XREFtrapz, *note dblquad: XREFdblquad, *note triplequad: XREFtriplequad. Sometimes one does not have the function, but only the raw (x, y) points from which to perform an integration. This can occur when collecting data in an experiment. The 'trapz' function can integrate these values as shown in the following example where "data" has been collected on the cosine function over the range [0, pi/2). x = 0:0.1:pi/2; # Uniformly spaced points y = cos (x); trapz (x, y) => 0.99666 The answer is reasonably close to the exact value of 1. Ordinary quadrature is sensitive to the characteristics of the integrand. Empirical integration depends not just on the integrand, but also on the particular points chosen to represent the function. Repeating the example above with the sine function over the range [0, pi/2) produces far inferior results. x = 0:0.1:pi/2; # Uniformly spaced points y = sin (x); trapz (x, y) => 0.92849 However, a slightly different choice of data points can change the result significantly. The same integration, with the same number of points, but spaced differently produces a more accurate answer. x = linspace (0, pi/2, 16); # Uniformly spaced, but including endpoint y = sin (x); trapz (x, y) => 0.99909 In general there may be no way of knowing the best distribution of points ahead of time. Or the points may come from an experiment where there is no freedom to select the best distribution. In any case, one must remain aware of this issue when using 'trapz'. -- Function File: Q = trapz (Y) -- Function File: Q = trapz (X, Y) -- Function File: Q = trapz (..., DIM) Numerically evaluate the integral of points Y using the trapezoidal method. 'trapz (Y)' computes the integral of Y along the first non-singleton dimension. When the argument X is omitted an equally spaced X vector with unit spacing (1) is assumed. 'trapz (X, Y)' evaluates the integral with respect to the spacing in X and the values in Y. This is useful if the points in Y have been sampled unevenly. If the optional DIM argument is given, operate along this dimension. Application Note: If X is not specified then unit spacing will be used. To scale the integral to the correct value you must multiply by the actual spacing value (deltaX). As an example, the integral of x^3 over the range [0, 1] is x^4/4 or 0.25. The following code uses 'trapz' to calculate the integral in three different ways. x = 0:0.1:1; y = x.^3; q = trapz (y) => q = 2.525 # No scaling q * 0.1 => q = 0.2525 # Approximation to integral by scaling trapz (x, y) => q = 0.2525 # Same result by specifying X See also: *note cumtrapz: XREFcumtrapz. -- Function File: Q = cumtrapz (Y) -- Function File: Q = cumtrapz (X, Y) -- Function File: Q = cumtrapz (..., DIM) Cumulative numerical integration of points Y using the trapezoidal method. 'cumtrapz (Y)' computes the cumulative integral of Y along the first non-singleton dimension. Where 'trapz' reports only the overall integral sum, 'cumtrapz' reports the current partial sum value at each point of Y. When the argument X is omitted an equally spaced X vector with unit spacing (1) is assumed. 'cumtrapz (X, Y)' evaluates the integral with respect to the spacing in X and the values in Y. This is useful if the points in Y have been sampled unevenly. If the optional DIM argument is given, operate along this dimension. Application Note: If X is not specified then unit spacing will be used. To scale the integral to the correct value you must multiply by the actual spacing value (deltaX). See also: *note trapz: XREFtrapz, *note cumsum: XREFcumsum.  File: octave.info, Node: Orthogonal Collocation, Next: Functions of Multiple Variables, Prev: Functions of One Variable, Up: Numerical Integration 23.2 Orthogonal Collocation =========================== -- Built-in Function: [R, AMAT, BMAT, Q] = colloc (N, "left", "right") Compute derivative and integral weight matrices for orthogonal collocation. Reference: J. Villadsen, M. L. Michelsen, 'Solution of Differential Equation Models by Polynomial Approximation'. Here is an example of using 'colloc' to generate weight matrices for solving the second order differential equation U' - ALPHA * U" = 0 with the boundary conditions U(0) = 0 and U(1) = 1. First, we can generate the weight matrices for N points (including the endpoints of the interval), and incorporate the boundary conditions in the right hand side (for a specific value of ALPHA). n = 7; alpha = 0.1; [r, a, b] = colloc (n-2, "left", "right"); at = a(2:n-1,2:n-1); bt = b(2:n-1,2:n-1); rhs = alpha * b(2:n-1,n) - a(2:n-1,n); Then the solution at the roots R is u = [ 0; (at - alpha * bt) \ rhs; 1] => [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]  File: octave.info, Node: Functions of Multiple Variables, Prev: Orthogonal Collocation, Up: Numerical Integration 23.3 Functions of Multiple Variables ==================================== Octave does not have built-in functions for computing the integral of functions of multiple variables directly. It is possible, however, to compute the integral of a function of multiple variables using the existing functions for one-dimensional integrals. To illustrate how the integration can be performed, we will integrate the function f(x, y) = sin(pi*x*y)*sqrt(x*y) for x and y between 0 and 1. The first approach creates a function that integrates f with respect to x, and then integrates that function with respect to y. Because 'quad' is written in Fortran it cannot be called recursively. This means that 'quad' cannot integrate a function that calls 'quad', and hence cannot be used to perform the double integration. Any of the other integrators, however, can be used which is what the following code demonstrates. function q = g(y) q = ones (size (y)); for i = 1:length (y) f = @(x) sin (pi*x.*y(i)) .* sqrt (x.*y(i)); q(i) = quadgk (f, 0, 1); endfor endfunction I = quadgk ("g", 0, 1) => 0.30022 The above process can be simplified with the 'dblquad' and 'triplequad' functions for integrals over two and three variables. For example: I = dblquad (@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1) => 0.30022 -- Function File: dblquad (F, XA, XB, YA, YB) -- Function File: dblquad (F, XA, XB, YA, YB, TOL) -- Function File: dblquad (F, XA, XB, YA, YB, TOL, QUADF) -- Function File: dblquad (F, XA, XB, YA, YB, TOL, QUADF, ...) Numerically evaluate the double integral of F. F is a function handle, inline function, or string containing the name of the function to evaluate. The function F must have the form z = f(x,y) where X is a vector and Y is a scalar. It should return a vector of the same length and orientation as X. XA, YA and XB, YB are the lower and upper limits of integration for x and y respectively. The underlying integrator determines whether infinite bounds are accepted. The optional argument TOL defines the absolute tolerance used to integrate each sub-integral. The default value is 1e^{-6}. The optional argument QUADF specifies which underlying integrator function to use. Any choice but 'quad' is available and the default is 'quadcc'. Additional arguments, are passed directly to F. To use the default value for TOL or QUADF one may pass ':' or an empty matrix ([]). See also: *note triplequad: XREFtriplequad, *note quad: XREFquad, *note quadv: XREFquadv, *note quadl: XREFquadl, *note quadgk: XREFquadgk, *note quadcc: XREFquadcc, *note trapz: XREFtrapz. -- Function File: triplequad (F, XA, XB, YA, YB, ZA, ZB) -- Function File: triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL) -- Function File: triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL, QUADF) -- Function File: triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL, QUADF, ...) Numerically evaluate the triple integral of F. F is a function handle, inline function, or string containing the name of the function to evaluate. The function F must have the form w = f(x,y,z) where either X or Y is a vector and the remaining inputs are scalars. It should return a vector of the same length and orientation as X or Y. XA, YA, ZA and XB, YB, ZB are the lower and upper limits of integration for x, y, and z respectively. The underlying integrator determines whether infinite bounds are accepted. The optional argument TOL defines the absolute tolerance used to integrate each sub-integral. The default value is 1e-6. The optional argument QUADF specifies which underlying integrator function to use. Any choice but 'quad' is available and the default is 'quadcc'. Additional arguments, are passed directly to F. To use the default value for TOL or QUADF one may pass ':' or an empty matrix ([]). See also: *note dblquad: XREFdblquad, *note quad: XREFquad, *note quadv: XREFquadv, *note quadl: XREFquadl, *note quadgk: XREFquadgk, *note quadcc: XREFquadcc, *note trapz: XREFtrapz. The above mentioned approach works, but is fairly slow, and that problem increases exponentially with the dimensionality of the integral. Another possible solution is to use Orthogonal Collocation as described in the previous section (*note Orthogonal Collocation::). The integral of a function f(x,y) for x and y between 0 and 1 can be approximated using n points by the sum over 'i=1:n' and 'j=1:n' of 'q(i)*q(j)*f(r(i),r(j))', where q and r is as returned by 'colloc (n)'. The generalization to more than two variables is straight forward. The following code computes the studied integral using n=8 points. f = @(x,y) sin (pi*x*y') .* sqrt (x*y'); n = 8; [t, ~, ~, q] = colloc (n); I = q'*f(t,t)*q; => 0.30022 It should be noted that the number of points determines the quality of the approximation. If the integration needs to be performed between a and b, instead of 0 and 1, then a change of variables is needed.  File: octave.info, Node: Differential Equations, Next: Optimization, Prev: Numerical Integration, Up: Top 24 Differential Equations ************************* Octave has built-in functions for solving ordinary differential equations, and differential-algebraic equations. All solvers are based on reliable ODE routines written in Fortran. * Menu: * Ordinary Differential Equations:: * Differential-Algebraic Equations::  File: octave.info, Node: Ordinary Differential Equations, Next: Differential-Algebraic Equations, Up: Differential Equations 24.1 Ordinary Differential Equations ==================================== The function 'lsode' can be used to solve ODEs of the form dx -- = f (x, t) dt using Hindmarsh's ODE solver LSODE. -- Built-in Function: [X, ISTATE, MSG] = lsode (FCN, X_0, T) -- Built-in Function: [X, ISTATE, MSG] = lsode (FCN, X_0, T, T_CRIT) Ordinary Differential Equation (ODE) solver. The set of differential equations to solve is dx -- = f (x, t) dt with x(t_0) = x_0 The solution is returned in the matrix X, with each row corresponding to an element of the vector T. The first element of T should be t_0 and should correspond to the initial state of the system X_0, so that the first row of the output is X_0. The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of right hand sides for the set of equations. The function must have the form XDOT = f (X, T) in which XDOT and X are vectors and T is a scalar. If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the Jacobian of f. The Jacobian function must have the form JAC = j (X, T) in which JAC is the matrix of partial derivatives | df_1 df_1 df_1 | | ---- ---- ... ---- | | dx_1 dx_2 dx_N | | | | df_2 df_2 df_2 | | ---- ---- ... ---- | df_i | dx_1 dx_2 dx_N | jac = ---- = | | dx_j | . . . . | | . . . . | | . . . . | | | | df_N df_N df_N | | ---- ---- ... ---- | | dx_1 dx_2 dx_N | The second and third arguments specify the initial state of the system, x_0, and the initial value of the independent variable t_0. The fourth argument is optional, and may be used to specify a set of times that the ODE solver should not integrate past. It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative. After a successful computation, the value of ISTATE will be 2 (consistent with the Fortran version of LSODE). If the computation is not successful, ISTATE will be something other than 2 and MSG will contain additional information. You can use the function 'lsode_options' to set optional parameters for 'lsode'. See also: *note daspk: XREFdaspk, *note dassl: XREFdassl, *note dasrt: XREFdasrt. -- Built-in Function: lsode_options () -- Built-in Function: val = lsode_options (OPT) -- Built-in Function: lsode_options (OPT, VAL) Query or set options for the function 'lsode'. When called with no arguments, the names of all available options and their current values are displayed. Given one argument, return the value of the option OPT. When called with two arguments, 'lsode_options' sets the option OPT to value VAL. Options include '"absolute tolerance"' Absolute tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector. '"relative tolerance"' Relative tolerance parameter. Unlike the absolute tolerance, this parameter may only be a scalar. The local error test applied at each integration step is abs (local error in x(i)) <= ... rtol * abs (y(i)) + atol(i) '"integration method"' A string specifying the method of integration to use to solve the ODE system. Valid values are "adams" "non-stiff" No Jacobian used (even if it is available). "bdf" "stiff" Use stiff backward differentiation formula (BDF) method. If a function to compute the Jacobian is not supplied, 'lsode' will compute a finite difference approximation of the Jacobian matrix. '"initial step size"' The step size to be attempted on the first step (default is determined automatically). '"maximum order"' Restrict the maximum order of the solution method. If using the Adams method, this option must be between 1 and 12. Otherwise, it must be between 1 and 5, inclusive. '"maximum step size"' Setting the maximum stepsize will avoid passing over very large regions (default is not specified). '"minimum step size"' The minimum absolute step size allowed (default is 0). '"step limit"' Maximum number of steps allowed (default is 100000). Here is an example of solving a set of three differential equations using 'lsode'. Given the function ## oregonator differential equation function xdot = f (x, t) xdot = zeros (3,1); xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) \ - 8.375e-06*x(1)^2); xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27; xdot(3) = 0.161*(x(1) - x(3)); endfunction and the initial condition 'x0 = [ 4; 1.1; 4 ]', the set of equations can be integrated using the command t = linspace (0, 500, 1000); y = lsode ("f", x0, t); If you try this, you will see that the value of the result changes dramatically between T = 0 and 5, and again around T = 305. A more efficient set of output points might be t = [0, logspace(-1, log10(303), 150), \ logspace(log10(304), log10(500), 150)]; See Alan C. Hindmarsh, 'ODEPACK, A Systematized Collection of ODE Solvers', in Scientific Computing, R. S. Stepleman, editor, (1983) for more information about the inner workings of 'lsode'. An m-file for the differential equation used above is included with the Octave distribution in the examples directory under the name 'oregonator.m'.  File: octave.info, Node: Differential-Algebraic Equations, Prev: Ordinary Differential Equations, Up: Differential Equations 24.2 Differential-Algebraic Equations ===================================== The function 'daspk' can be used to solve DAEs of the form 0 = f (x-dot, x, t), x(t=0) = x_0, x-dot(t=0) = x-dot_0 where x-dot is the derivative of x. The equation is solved using Petzold's DAE solver DASPK. -- Built-in Function: [X, XDOT, ISTATE, MSG] = daspk (FCN, X_0, XDOT_0, T, T_CRIT) Solve the set of differential-algebraic equations 0 = f (x, xdot, t) with x(t_0) = x_0, xdot(t_0) = xdot_0 The solution is returned in the matrices X and XDOT, with each row in the result matrices corresponding to one of the elements in the vector T. The first element of T should be t_0 and correspond to the initial state of the system X_0 and its derivative XDOT_0, so that the first row of the output X is X_0 and the first row of the output XDOT is XDOT_0. The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of residuals for the set of equations. It must have the form RES = f (X, XDOT, T) in which X, XDOT, and RES are vectors, and T is a scalar. If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the modified Jacobian df df jac = -- + c ------ dx d xdot The modified Jacobian function must have the form JAC = j (X, XDOT, T, C) The second and third arguments to 'daspk' specify the initial condition of the states and their derivatives, and the fourth argument specifies a vector of output times at which the solution is desired, including the time corresponding to the initial condition. The set of initial states and derivatives are not strictly required to be consistent. If they are not consistent, you must use the 'daspk_options' function to provide additional information so that 'daspk' can compute a consistent starting point. The fifth argument is optional, and may be used to specify a set of times that the DAE solver should not integrate past. It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative. After a successful computation, the value of ISTATE will be greater than zero (consistent with the Fortran version of DASPK). If the computation is not successful, the value of ISTATE will be less than zero and MSG will contain additional information. You can use the function 'daspk_options' to set optional parameters for 'daspk'. See also: *note dassl: XREFdassl. -- Built-in Function: daspk_options () -- Built-in Function: val = daspk_options (OPT) -- Built-in Function: daspk_options (OPT, VAL) Query or set options for the function 'daspk'. When called with no arguments, the names of all available options and their current values are displayed. Given one argument, return the value of the option OPT. When called with two arguments, 'daspk_options' sets the option OPT to value VAL. Options include '"absolute tolerance"' Absolute tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector, and the relative tolerance must also be a vector of the same length. '"relative tolerance"' Relative tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector, and the absolute tolerance must also be a vector of the same length. The local error test applied at each integration step is abs (local error in x(i)) <= rtol(i) * abs (Y(i)) + atol(i) '"compute consistent initial condition"' Denoting the differential variables in the state vector by 'Y_d' and the algebraic variables by 'Y_a', 'ddaspk' can solve one of two initialization problems: 1. Given Y_d, calculate Y_a and Y'_d 2. Given Y', calculate Y. In either case, initial values for the given components are input, and initial guesses for the unknown components must also be provided as input. Set this option to 1 to solve the first problem, or 2 to solve the second (the default is 0, so you must provide a set of initial conditions that are consistent). If this option is set to a nonzero value, you must also set the "algebraic variables" option to declare which variables in the problem are algebraic. '"use initial condition heuristics"' Set to a nonzero value to use the initial condition heuristics options described below. '"initial condition heuristics"' A vector of the following parameters that can be used to control the initial condition calculation. 'MXNIT' Maximum number of Newton iterations (default is 5). 'MXNJ' Maximum number of Jacobian evaluations (default is 6). 'MXNH' Maximum number of values of the artificial stepsize parameter to be tried if the "compute consistent initial condition" option has been set to 1 (default is 5). Note that the maximum total number of Newton iterations allowed is 'MXNIT*MXNJ*MXNH' if the "compute consistent initial condition" option has been set to 1 and 'MXNIT*MXNJ' if it is set to 2. 'LSOFF' Set to a nonzero value to disable the linesearch algorithm (default is 0). 'STPTOL' Minimum scaled step in linesearch algorithm (default is eps^(2/3)). 'EPINIT' Swing factor in the Newton iteration convergence test. The test is applied to the residual vector, premultiplied by the approximate Jacobian. For convergence, the weighted RMS norm of this vector (scaled by the error weights) must be less than 'EPINIT*EPCON', where 'EPCON' = 0.33 is the analogous test constant used in the time steps. The default is 'EPINIT' = 0.01. '"print initial condition info"' Set this option to a nonzero value to display detailed information about the initial condition calculation (default is 0). '"exclude algebraic variables from error test"' Set to a nonzero value to exclude algebraic variables from the error test. You must also set the "algebraic variables" option to declare which variables in the problem are algebraic (default is 0). '"algebraic variables"' A vector of the same length as the state vector. A nonzero element indicates that the corresponding element of the state vector is an algebraic variable (i.e., its derivative does not appear explicitly in the equation set). This option is required by the "compute consistent initial condition" and "exclude algebraic variables from error test" options. '"enforce inequality constraints"' Set to one of the following values to enforce the inequality constraints specified by the "inequality constraint types" option (default is 0). 1. To have constraint checking only in the initial condition calculation. 2. To enforce constraint checking during the integration. 3. To enforce both options 1 and 2. '"inequality constraint types"' A vector of the same length as the state specifying the type of inequality constraint. Each element of the vector corresponds to an element of the state and should be assigned one of the following codes -2 Less than zero. -1 Less than or equal to zero. 0 Not constrained. 1 Greater than or equal to zero. 2 Greater than zero. This option only has an effect if the "enforce inequality constraints" option is nonzero. '"initial step size"' Differential-algebraic problems may occasionally suffer from severe scaling difficulties on the first step. If you know a great deal about the scaling of your problem, you can help to alleviate this problem by specifying an initial stepsize (default is computed automatically). '"maximum order"' Restrict the maximum order of the solution method. This option must be between 1 and 5, inclusive (default is 5). '"maximum step size"' Setting the maximum stepsize will avoid passing over very large regions (default is not specified). Octave also includes DASSL, an earlier version of DASPK, and DASRT, which can be used to solve DAEs with constraints (stopping conditions). -- Built-in Function: [X, XDOT, ISTATE, MSG] = dassl (FCN, X_0, XDOT_0, T, T_CRIT) Solve the set of differential-algebraic equations 0 = f (x, xdot, t) with x(t_0) = x_0, xdot(t_0) = xdot_0 The solution is returned in the matrices X and XDOT, with each row in the result matrices corresponding to one of the elements in the vector T. The first element of T should be t_0 and correspond to the initial state of the system X_0 and its derivative XDOT_0, so that the first row of the output X is X_0 and the first row of the output XDOT is XDOT_0. The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of residuals for the set of equations. It must have the form RES = f (X, XDOT, T) in which X, XDOT, and RES are vectors, and T is a scalar. If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the modified Jacobian df df jac = -- + c ------ dx d xdot The modified Jacobian function must have the form JAC = j (X, XDOT, T, C) The second and third arguments to 'dassl' specify the initial condition of the states and their derivatives, and the fourth argument specifies a vector of output times at which the solution is desired, including the time corresponding to the initial condition. The set of initial states and derivatives are not strictly required to be consistent. In practice, however, DASSL is not very good at determining a consistent set for you, so it is best if you ensure that the initial values result in the function evaluating to zero. The fifth argument is optional, and may be used to specify a set of times that the DAE solver should not integrate past. It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative. After a successful computation, the value of ISTATE will be greater than zero (consistent with the Fortran version of DASSL). If the computation is not successful, the value of ISTATE will be less than zero and MSG will contain additional information. You can use the function 'dassl_options' to set optional parameters for 'dassl'. See also: *note daspk: XREFdaspk, *note dasrt: XREFdasrt, *note lsode: XREFlsode. -- Built-in Function: dassl_options () -- Built-in Function: val = dassl_options (OPT) -- Built-in Function: dassl_options (OPT, VAL) Query or set options for the function 'dassl'. When called with no arguments, the names of all available options and their current values are displayed. Given one argument, return the value of the option OPT. When called with two arguments, 'dassl_options' sets the option OPT to value VAL. Options include '"absolute tolerance"' Absolute tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector, and the relative tolerance must also be a vector of the same length. '"relative tolerance"' Relative tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector, and the absolute tolerance must also be a vector of the same length. The local error test applied at each integration step is abs (local error in x(i)) <= rtol(i) * abs (Y(i)) + atol(i) '"compute consistent initial condition"' If nonzero, 'dassl' will attempt to compute a consistent set of initial conditions. This is generally not reliable, so it is best to provide a consistent set and leave this option set to zero. '"enforce nonnegativity constraints"' If you know that the solutions to your equations will always be non-negative, it may help to set this parameter to a nonzero value. However, it is probably best to try leaving this option set to zero first, and only setting it to a nonzero value if that doesn't work very well. '"initial step size"' Differential-algebraic problems may occasionally suffer from severe scaling difficulties on the first step. If you know a great deal about the scaling of your problem, you can help to alleviate this problem by specifying an initial stepsize. '"maximum order"' Restrict the maximum order of the solution method. This option must be between 1 and 5, inclusive. '"maximum step size"' Setting the maximum stepsize will avoid passing over very large regions (default is not specified). '"step limit"' Maximum number of integration steps to attempt on a single call to the underlying Fortran code. -- Built-in Function: [X, XDOT, T_OUT, ISTAT, MSG] = dasrt (FCN, [], X_0, XDOT_0, T) -- Built-in Function: ... = dasrt (FCN, G, X_0, XDOT_0, T) -- Built-in Function: ... = dasrt (FCN, [], X_0, XDOT_0, T, T_CRIT) -- Built-in Function: ... = dasrt (FCN, G, X_0, XDOT_0, T, T_CRIT) Solve the set of differential-algebraic equations 0 = f (x, xdot, t) with x(t_0) = x_0, xdot(t_0) = xdot_0 with functional stopping criteria (root solving). The solution is returned in the matrices X and XDOT, with each row in the result matrices corresponding to one of the elements in the vector T_OUT. The first element of T should be t_0 and correspond to the initial state of the system X_0 and its derivative XDOT_0, so that the first row of the output X is X_0 and the first row of the output XDOT is XDOT_0. The vector T provides an upper limit on the length of the integration. If the stopping condition is met, the vector T_OUT will be shorter than T, and the final element of T_OUT will be the point at which the stopping condition was met, and may not correspond to any element of the vector T. The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of residuals for the set of equations. It must have the form RES = f (X, XDOT, T) in which X, XDOT, and RES are vectors, and T is a scalar. If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the modified Jacobian df df jac = -- + c ------ dx d xdot The modified Jacobian function must have the form JAC = j (X, XDOT, T, C) The optional second argument names a function that defines the constraint functions whose roots are desired during the integration. This function must have the form G_OUT = g (X, T) and return a vector of the constraint function values. If the value of any of the constraint functions changes sign, DASRT will attempt to stop the integration at the point of the sign change. If the name of the constraint function is omitted, 'dasrt' solves the same problem as 'daspk' or 'dassl'. Note that because of numerical errors in the constraint functions due to round-off and integration error, DASRT may return false roots, or return the same root at two or more nearly equal values of T. If such false roots are suspected, the user should consider smaller error tolerances or higher precision in the evaluation of the constraint functions. If a root of some constraint function defines the end of the problem, the input to DASRT should nevertheless allow integration to a point slightly past that root, so that DASRT can locate the root by interpolation. The third and fourth arguments to 'dasrt' specify the initial condition of the states and their derivatives, and the fourth argument specifies a vector of output times at which the solution is desired, including the time corresponding to the initial condition. The set of initial states and derivatives are not strictly required to be consistent. In practice, however, DASSL is not very good at determining a consistent set for you, so it is best if you ensure that the initial values result in the function evaluating to zero. The sixth argument is optional, and may be used to specify a set of times that the DAE solver should not integrate past. It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative. After a successful computation, the value of ISTATE will be greater than zero (consistent with the Fortran version of DASSL). If the computation is not successful, the value of ISTATE will be less than zero and MSG will contain additional information. You can use the function 'dasrt_options' to set optional parameters for 'dasrt'. See also: *note dasrt_options: XREFdasrt_options, *note daspk: XREFdaspk, *note dasrt: XREFdasrt, *note lsode: XREFlsode. -- Built-in Function: dasrt_options () -- Built-in Function: val = dasrt_options (OPT) -- Built-in Function: dasrt_options (OPT, VAL) Query or set options for the function 'dasrt'. When called with no arguments, the names of all available options and their current values are displayed. Given one argument, return the value of the option OPT. When called with two arguments, 'dasrt_options' sets the option OPT to value VAL. Options include '"absolute tolerance"' Absolute tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector, and the relative tolerance must also be a vector of the same length. '"relative tolerance"' Relative tolerance. May be either vector or scalar. If a vector, it must match the dimension of the state vector, and the absolute tolerance must also be a vector of the same length. The local error test applied at each integration step is abs (local error in x(i)) <= ... rtol(i) * abs (Y(i)) + atol(i) '"initial step size"' Differential-algebraic problems may occasionally suffer from severe scaling difficulties on the first step. If you know a great deal about the scaling of your problem, you can help to alleviate this problem by specifying an initial stepsize. '"maximum order"' Restrict the maximum order of the solution method. This option must be between 1 and 5, inclusive. '"maximum step size"' Setting the maximum stepsize will avoid passing over very large regions. '"step limit"' Maximum number of integration steps to attempt on a single call to the underlying Fortran code. See K. E. Brenan, et al., 'Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations', North-Holland (1989) for more information about the implementation of DASSL.  File: octave.info, Node: Optimization, Next: Statistics, Prev: Differential Equations, Up: Top 25 Optimization *************** Octave comes with support for solving various kinds of optimization problems. Specifically Octave can solve problems in Linear Programming, Quadratic Programming, Nonlinear Programming, and Linear Least Squares Minimization. * Menu: * Linear Programming:: * Quadratic Programming:: * Nonlinear Programming:: * Linear Least Squares::  File: octave.info, Node: Linear Programming, Next: Quadratic Programming, Up: Optimization 25.1 Linear Programming ======================= Octave can solve Linear Programming problems using the 'glpk' function. That is, Octave can solve min C'*x subject to the linear constraints A*x = b where x >= 0. The 'glpk' function also supports variations of this problem. -- Function File: [XOPT, FMIN, ERRNUM, EXTRA] = glpk (C, A, B, LB, UB, CTYPE, VARTYPE, SENSE, PARAM) Solve a linear program using the GNU GLPK library. Given three arguments, 'glpk' solves the following standard LP: min C'*x subject to A*x = b x >= 0 but may also solve problems of the form [ min | max ] C'*x subject to A*x [ "=" | "<=" | ">=" ] b x >= LB x <= UB Input arguments: C A column array containing the objective function coefficients. A A matrix containing the constraints coefficients. B A column array containing the right-hand side value for each constraint in the constraint matrix. LB An array containing the lower bound on each of the variables. If LB is not supplied, the default lower bound for the variables is zero. UB An array containing the upper bound on each of the variables. If UB is not supplied, the default upper bound is assumed to be infinite. CTYPE An array of characters containing the sense of each constraint in the constraint matrix. Each element of the array may be one of the following values "F" A free (unbounded) constraint (the constraint is ignored). "U" An inequality constraint with an upper bound ('A(i,:)*x <= b(i)'). "S" An equality constraint ('A(i,:)*x = b(i)'). "L" An inequality with a lower bound ('A(i,:)*x >= b(i)'). "D" An inequality constraint with both upper and lower bounds ('A(i,:)*x >= -b(i)') _and_ ('A(i,:)*x <= b(i)'). VARTYPE A column array containing the types of the variables. "C" A continuous variable. "I" An integer variable. SENSE If SENSE is 1, the problem is a minimization. If SENSE is -1, the problem is a maximization. The default value is 1. PARAM A structure containing the following parameters used to define the behavior of solver. Missing elements in the structure take on default values, so you only need to set the elements that you wish to change from the default. Integer parameters: 'msglev (default: 1)' Level of messages output by solver routines: 0 ('GLP_MSG_OFF') No output. 1 ('GLP_MSG_ERR') Error and warning messages only. 2 ('GLP_MSG_ON') Normal output. 3 ('GLP_MSG_ALL') Full output (includes informational messages). 'scale (default: 16)' Scaling option. The values can be combined with the bitwise OR operator and may be the following: 1 ('GLP_SF_GM') Geometric mean scaling. 16 ('GLP_SF_EQ') Equilibration scaling. 32 ('GLP_SF_2N') Round scale factors to power of two. 64 ('GLP_SF_SKIP') Skip if problem is well scaled. Alternatively, a value of 128 ('GLP_SF_AUTO') may be also specified, in which case the routine chooses the scaling options automatically. 'dual (default: 1)' Simplex method option: 1 ('GLP_PRIMAL') Use two-phase primal simplex. 2 ('GLP_DUALP') Use two-phase dual simplex, and if it fails, switch to the primal simplex. 3 ('GLP_DUAL') Use two-phase dual simplex. 'price (default: 34)' Pricing option (for both primal and dual simplex): 17 ('GLP_PT_STD') Textbook pricing. 34 ('GLP_PT_PSE') Steepest edge pricing. 'itlim (default: intmax)' Simplex iterations limit. It is decreased by one each time when one simplex iteration has been performed, and reaching zero value signals the solver to stop the search. 'outfrq (default: 200)' Output frequency, in iterations. This parameter specifies how frequently the solver sends information about the solution to the standard output. 'branch (default: 4)' Branching technique option (for MIP only): 1 ('GLP_BR_FFV') First fractional variable. 2 ('GLP_BR_LFV') Last fractional variable. 3 ('GLP_BR_MFV') Most fractional variable. 4 ('GLP_BR_DTH') Heuristic by Driebeck and Tomlin. 5 ('GLP_BR_PCH') Hybrid pseudocost heuristic. 'btrack (default: 4)' Backtracking technique option (for MIP only): 1 ('GLP_BT_DFS') Depth first search. 2 ('GLP_BT_BFS') Breadth first search. 3 ('GLP_BT_BLB') Best local bound. 4 ('GLP_BT_BPH') Best projection heuristic. 'presol (default: 1)' If this flag is set, the simplex solver uses the built-in LP presolver. Otherwise the LP presolver is not used. 'lpsolver (default: 1)' Select which solver to use. If the problem is a MIP problem this flag will be ignored. 1 Revised simplex method. 2 Interior point method. 'rtest (default: 34)' Ratio test technique: 17 ('GLP_RT_STD') Standard ("textbook"). 34 ('GLP_RT_HAR') Harris' two-pass ratio test. 'tmlim (default: intmax)' Searching time limit, in milliseconds. 'outdly (default: 0)' Output delay, in seconds. This parameter specifies how long the solver should delay sending information about the solution to the standard output. 'save (default: 0)' If this parameter is nonzero, save a copy of the problem in CPLEX LP format to the file '"outpb.lp"'. There is currently no way to change the name of the output file. Real parameters: 'tolbnd (default: 1e-7)' Relative tolerance used to check if the current basic solution is primal feasible. It is not recommended that you change this parameter unless you have a detailed understanding of its purpose. 'toldj (default: 1e-7)' Absolute tolerance used to check if the current basic solution is dual feasible. It is not recommended that you change this parameter unless you have a detailed understanding of its purpose. 'tolpiv (default: 1e-10)' Relative tolerance used to choose eligible pivotal elements of the simplex table. It is not recommended that you change this parameter unless you have a detailed understanding of its purpose. 'objll (default: -DBL_MAX)' Lower limit of the objective function. If the objective function reaches this limit and continues decreasing, the solver stops the search. This parameter is used in the dual simplex method only. 'objul (default: +DBL_MAX)' Upper limit of the objective function. If the objective function reaches this limit and continues increasing, the solver stops the search. This parameter is used in the dual simplex only. 'tolint (default: 1e-5)' Relative tolerance used to check if the current basic solution is integer feasible. It is not recommended that you change this parameter unless you have a detailed understanding of its purpose. 'tolobj (default: 1e-7)' Relative tolerance used to check if the value of the objective function is not better than in the best known integer feasible solution. It is not recommended that you change this parameter unless you have a detailed understanding of its purpose. Output values: XOPT The optimizer (the value of the decision variables at the optimum). FOPT The optimum value of the objective function. ERRNUM Error code. 0 No error. 1 ('GLP_EBADB') Invalid basis. 2 ('GLP_ESING') Singular matrix. 3 ('GLP_ECOND') Ill-conditioned matrix. 4 ('GLP_EBOUND') Invalid bounds. 5 ('GLP_EFAIL') Solver failed. 6 ('GLP_EOBJLL') Objective function lower limit reached. 7 ('GLP_EOBJUL') Objective function upper limit reached. 8 ('GLP_EITLIM') Iterations limit exhausted. 9 ('GLP_ETMLIM') Time limit exhausted. 10 ('GLP_ENOPFS') No primal feasible solution. 11 ('GLP_ENODFS') No dual feasible solution. 12 ('GLP_EROOT') Root LP optimum not provided. 13 ('GLP_ESTOP') Search terminated by application. 14 ('GLP_EMIPGAP') Relative MIP gap tolerance reached. 15 ('GLP_ENOFEAS') No primal/dual feasible solution. 16 ('GLP_ENOCVG') No convergence. 17 ('GLP_EINSTAB') Numerical instability. 18 ('GLP_EDATA') Invalid data. 19 ('GLP_ERANGE') Result out of range. EXTRA A data structure containing the following fields: 'lambda' Dual variables. 'redcosts' Reduced Costs. 'time' Time (in seconds) used for solving LP/MIP problem. 'status' Status of the optimization. 1 ('GLP_UNDEF') Solution status is undefined. 2 ('GLP_FEAS') Solution is feasible. 3 ('GLP_INFEAS') Solution is infeasible. 4 ('GLP_NOFEAS') Problem has no feasible solution. 5 ('GLP_OPT') Solution is optimal. 6 ('GLP_UNBND') Problem has no unbounded solution. Example: c = [10, 6, 4]'; A = [ 1, 1, 1; 10, 4, 5; 2, 2, 6]; b = [100, 600, 300]'; lb = [0, 0, 0]'; ub = []; ctype = "UUU"; vartype = "CCC"; s = -1; param.msglev = 1; param.itlim = 100; [xmin, fmin, status, extra] = ... glpk (c, A, b, lb, ub, ctype, vartype, s, param);  File: octave.info, Node: Quadratic Programming, Next: Nonlinear Programming, Prev: Linear Programming, Up: Optimization 25.2 Quadratic Programming ========================== Octave can also solve Quadratic Programming problems, this is min 0.5 x'*H*x + x'*q subject to A*x = b lb <= x <= ub A_lb <= A_in*x <= A_ub -- Function File: [X, OBJ, INFO, LAMBDA] = qp (X0, H) -- Function File: [X, OBJ, INFO, LAMBDA] = qp (X0, H, Q) -- Function File: [X, OBJ, INFO, LAMBDA] = qp (X0, H, Q, A, B) -- Function File: [X, OBJ, INFO, LAMBDA] = qp (X0, H, Q, A, B, LB, UB) -- Function File: [X, OBJ, INFO, LAMBDA] = qp (X0, H, Q, A, B, LB, UB, A_LB, A_IN, A_UB) -- Function File: [X, OBJ, INFO, LAMBDA] = qp (..., OPTIONS) Solve a quadratic program (QP). Solve the quadratic program defined by min 0.5 x'*H*x + x'*q x subject to A*x = b lb <= x <= ub A_lb <= A_in*x <= A_ub using a null-space active-set method. Any bound (A, B, LB, UB, A_LB, A_UB) may be set to the empty matrix ('[]') if not present. If the initial guess is feasible the algorithm is faster. OPTIONS An optional structure containing the following parameter(s) used to define the behavior of the solver. Missing elements in the structure take on default values, so you only need to set the elements that you wish to change from the default. 'MaxIter (default: 200)' Maximum number of iterations. INFO Structure containing run-time information about the algorithm. The following fields are defined: 'solveiter' The number of iterations required to find the solution. 'info' An integer indicating the status of the solution. 0 The problem is feasible and convex. Global solution found. 1 The problem is not convex. Local solution found. 2 The problem is not convex and unbounded. 3 Maximum number of iterations reached. 6 The problem is infeasible. -- Function File: X = pqpnonneg (C, D) -- Function File: X = pqpnonneg (C, D, X0) -- Function File: [X, MINVAL] = pqpnonneg (...) -- Function File: [X, MINVAL, EXITFLAG] = pqpnonneg (...) -- Function File: [X, MINVAL, EXITFLAG, OUTPUT] = pqpnonneg (...) -- Function File: [X, MINVAL, EXITFLAG, OUTPUT, LAMBDA] = pqpnonneg (...) Minimize '1/2*x'*c*x + d'*x' subject to 'X >= 0'. C ## and D must be real, and C must be symmetric and positive definite. X0 is an optional initial guess for X. Outputs: * minval The minimum attained model value, 1/2*xmin'*c*xmin + d'*xmin * exitflag An indicator of convergence. 0 indicates that the iteration count was exceeded, and therefore convergence was not reached; >0 indicates that the algorithm converged. (The algorithm is stable and will converge given enough iterations.) * output A structure with two fields: * "algorithm": The algorithm used ("nnls") * "iterations": The number of iterations taken. * lambda Not implemented. See also: *note optimset: XREFoptimset, *note lsqnonneg: XREFlsqnonneg, *note qp: XREFqp.  File: octave.info, Node: Nonlinear Programming, Next: Linear Least Squares, Prev: Quadratic Programming, Up: Optimization 25.3 Nonlinear Programming ========================== Octave can also perform general nonlinear minimization using a successive quadratic programming solver. -- Function File: [X, OBJ, INFO, ITER, NF, LAMBDA] = sqp (X0, PHI) -- Function File: [...] = sqp (X0, PHI, G) -- Function File: [...] = sqp (X0, PHI, G, H) -- Function File: [...] = sqp (X0, PHI, G, H, LB, UB) -- Function File: [...] = sqp (X0, PHI, G, H, LB, UB, MAXITER) -- Function File: [...] = sqp (X0, PHI, G, H, LB, UB, MAXITER, TOL) Minimize an objective function using sequential quadratic programming (SQP). Solve the nonlinear program min phi (x) x subject to g(x) = 0 h(x) >= 0 lb <= x <= ub using a sequential quadratic programming method. The first argument is the initial guess for the vector X0. The second argument is a function handle pointing to the objective function PHI. The objective function must accept one vector argument and return a scalar. The second argument may also be a 2- or 3-element cell array of function handles. The first element should point to the objective function, the second should point to a function that computes the gradient of the objective function, and the third should point to a function that computes the Hessian of the objective function. If the gradient function is not supplied, the gradient is computed by finite differences. If the Hessian function is not supplied, a BFGS update formula is used to approximate the Hessian. When supplied, the gradient function 'PHI{2}' must accept one vector argument and return a vector. When supplied, the Hessian function 'PHI{3}' must accept one vector argument and return a matrix. The third and fourth arguments G and H are function handles pointing to functions that compute the equality constraints and the inequality constraints, respectively. If the problem does not have equality (or inequality) constraints, then use an empty matrix ([]) for G (or H). When supplied, these equality and inequality constraint functions must accept one vector argument and return a vector. The third and fourth arguments may also be 2-element cell arrays of function handles. The first element should point to the constraint function and the second should point to a function that computes the gradient of the constraint function: [ d f(x) d f(x) d f(x) ] transpose ( [ ------ ----- ... ------ ] ) [ dx_1 dx_2 dx_N ] The fifth and sixth arguments, LB and UB, contain lower and upper bounds on X. These must be consistent with the equality and inequality constraints G and H. If the arguments are vectors then X(i) is bound by LB(i) and UB(i). A bound can also be a scalar in which case all elements of X will share the same bound. If only one bound (lb, ub) is specified then the other will default to (-REALMAX, +REALMAX). The seventh argument MAXITER specifies the maximum number of iterations. The default value is 100. The eighth argument TOL specifies the tolerance for the stopping criteria. The default value is 'sqrt (eps)'. The value returned in INFO may be one of the following: 101 The algorithm terminated normally. All constraints meet the specified tolerance. 102 The BFGS update failed. 103 The maximum number of iterations was reached. 104 The stepsize has become too small, i.e., delta X, is less than 'TOL * norm (x)'. An example of calling 'sqp': function r = g (x) r = [ sumsq(x)-10; x(2)*x(3)-5*x(4)*x(5); x(1)^3+x(2)^3+1 ]; endfunction function obj = phi (x) obj = exp (prod (x)) - 0.5*(x(1)^3+x(2)^3+1)^2; endfunction x0 = [-1.8; 1.7; 1.9; -0.8; -0.8]; [x, obj, info, iter, nf, lambda] = sqp (x0, @phi, @g, []) x = -1.71714 1.59571 1.82725 -0.76364 -0.76364 obj = 0.053950 info = 101 iter = 8 nf = 10 lambda = -0.0401627 0.0379578 -0.0052227 See also: *note qp: XREFqp.  File: octave.info, Node: Linear Least Squares, Prev: Nonlinear Programming, Up: Optimization 25.4 Linear Least Squares ========================= Octave also supports linear least squares minimization. That is, Octave can find the parameter b such that the model y = x*b fits data (x,y) as well as possible, assuming zero-mean Gaussian noise. If the noise is assumed to be isotropic the problem can be solved using the '\' or '/' operators, or the 'ols' function. In the general case where the noise is assumed to be anisotropic the 'gls' is needed. -- Function File: [BETA, SIGMA, R] = ols (Y, X) Ordinary least squares estimation. OLS applies to the multivariate model y = x*b + e with mean (e) = 0 and cov (vec (e)) = kron (s, I). where y is a t by p matrix, x is a t by k matrix, b is a k by p matrix, and e is a t by p matrix. Each row of Y and X is an observation and each column a variable. The return values BETA, SIGMA, and R are defined as follows. BETA The OLS estimator for b. BETA is calculated directly via 'inv (x'*x) * x' * y' if the matrix 'x'*x' is of full rank. Otherwise, 'BETA = pinv (X) * Y' where 'pinv (X)' denotes the pseudoinverse of X. SIGMA The OLS estimator for the matrix S, SIGMA = (Y-X*BETA)' * (Y-X*BETA) / (T-rank(X)) R The matrix of OLS residuals, 'R = Y - X*BETA'. See also: *note gls: XREFgls, *note pinv: XREFpinv. -- Function File: [BETA, V, R] = gls (Y, X, O) Generalized least squares model. Perform a generalized least squares estimation for the multivariate model y = x*b + e with mean (e) = 0 and cov (vec (e)) = (s^2) o, where y is a t by p matrix, x is a t by k matrix, b is a k by p matrix, e is a t by p matrix, and o is a t*p by t*p matrix. Each row of Y and X is an observation and each column a variable. The return values BETA, V, and R are defined as follows. BETA The GLS estimator for b. V The GLS estimator for s^2. R The matrix of GLS residuals, r = y - x*beta. See also: *note ols: XREFols. -- Function File: X = lsqnonneg (C, D) -- Function File: X = lsqnonneg (C, D, X0) -- Function File: X = lsqnonneg (C, D, X0, OPTIONS) -- Function File: [X, RESNORM] = lsqnonneg (...) -- Function File: [X, RESNORM, RESIDUAL] = lsqnonneg (...) -- Function File: [X, RESNORM, RESIDUAL, EXITFLAG] = lsqnonneg (...) -- Function File: [X, RESNORM, RESIDUAL, EXITFLAG, OUTPUT] = lsqnonneg (...) -- Function File: [X, RESNORM, RESIDUAL, EXITFLAG, OUTPUT, LAMBDA] = lsqnonneg (...) Minimize 'norm (C*X - d)' subject to 'X >= 0'. C and D must be real. X0 is an optional initial guess for X. Currently, 'lsqnonneg' recognizes these options: "MaxIter", "TolX". For a description of these options, see *note optimset: XREFoptimset. Outputs: * resnorm The squared 2-norm of the residual: norm (C*X-D)^2 * residual The residual: D-C*X * exitflag An indicator of convergence. 0 indicates that the iteration count was exceeded, and therefore convergence was not reached; >0 indicates that the algorithm converged. (The algorithm is stable and will converge given enough iterations.) * output A structure with two fields: * "algorithm": The algorithm used ("nnls") * "iterations": The number of iterations taken. * lambda Not implemented. See also: *note optimset: XREFoptimset, *note pqpnonneg: XREFpqpnonneg, *note lscov: XREFlscov. -- Function File: X = lscov (A, B) -- Function File: X = lscov (A, B, V) -- Function File: X = lscov (A, B, V, ALG) -- Function File: [X, STDX, MSE, S] = lscov (...) Compute a generalized linear least squares fit. Estimate X under the model B = AX + W, where the noise W is assumed to follow a normal distribution with covariance matrix {\sigma^2} V. If the size of the coefficient matrix A is n-by-p, the size of the vector/array of constant terms B must be n-by-k. The optional input argument V may be a n-by-1 vector of positive weights (inverse variances), or a n-by-n symmetric positive semidefinite matrix representing the covariance of B. If V is not supplied, the ordinary least squares solution is returned. The ALG input argument, a guidance on solution method to use, is currently ignored. Besides the least-squares estimate matrix X (p-by-k), the function also returns STDX (p-by-k), the error standard deviation of estimated X; MSE (k-by-1), the estimated data error covariance scale factors (\sigma^2); and S (p-by-p, or p-by-p-by-k if k > 1), the error covariance of X. Reference: Golub and Van Loan (1996), 'Matrix Computations (3rd Ed.)', Johns Hopkins, Section 5.6.3 See also: *note ols: XREFols, *note gls: XREFgls, *note lsqnonneg: XREFlsqnonneg. -- Function File: optimset () -- Function File: OPTIONS = optimset () -- Function File: OPTIONS = optimset (PAR, VAL, ...) -- Function File: OPTIONS = optimset (OLD, PAR, VAL, ...) -- Function File: OPTIONS = optimset (OLD, NEW) Create options structure for optimization functions. When called without any input or output arguments, 'optimset' prints a list of all valid optimization parameters. When called with one output and no inputs, return an options structure with all valid option parameters initialized to '[]'. When called with a list of parameter/value pairs, return an options structure with only the named parameters initialized. When the first input is an existing options structure OLD, the values are updated from either the PAR/VAL list or from the options structure NEW. Valid parameters are: AutoScaling ComplexEqn Display Request verbose display of results from optimizations. Values are: "off" [default] No display. "iter" Display intermediate results for every loop iteration. "final" Display the result of the final loop iteration. "notify" Display the result of the final loop iteration if the function has failed to converge. FinDiffType FunValCheck When enabled, display an error if the objective function returns an invalid value (a complex number, NaN, or Inf). Must be set to "on" or "off" [default]. Note: the functions 'fzero' and 'fminbnd' correctly handle Inf values and only complex values or NaN will cause an error in this case. GradObj When set to "on", the function to be minimized must return a second argument which is the gradient, or first derivative, of the function at the point X. If set to "off" [default], the gradient is computed via finite differences. Jacobian When set to "on", the function to be minimized must return a second argument which is the Jacobian, or first derivative, of the function at the point X. If set to "off" [default], the Jacobian is computed via finite differences. MaxFunEvals Maximum number of function evaluations before optimization stops. Must be a positive integer. MaxIter Maximum number of algorithm iterations before optimization stops. Must be a positive integer. OutputFcn A user-defined function executed once per algorithm iteration. TolFun Termination criterion for the function output. If the difference in the calculated objective function between one algorithm iteration and the next is less than 'TolFun' the optimization stops. Must be a positive scalar. TolX Termination criterion for the function input. If the difference in X, the current search point, between one algorithm iteration and the next is less than 'TolX' the optimization stops. Must be a positive scalar. TypicalX Updating See also: *note optimget: XREFoptimget. -- Function File: optimget (OPTIONS, PARNAME) -- Function File: optimget (OPTIONS, PARNAME, DEFAULT) Return the specific option PARNAME from the optimization options structure OPTIONS created by 'optimset'. If PARNAME is not defined then return DEFAULT if supplied, otherwise return an empty matrix. See also: *note optimset: XREFoptimset.  File: octave.info, Node: Statistics, Next: Sets, Prev: Optimization, Up: Top 26 Statistics ************* Octave has support for various statistical methods. This includes basic descriptive statistics, probability distributions, statistical tests, random number generation, and much more. The functions that analyze data all assume that multi-dimensional data is arranged in a matrix where each row is an observation, and each column is a variable. Thus, the matrix defined by a = [ 0.9, 0.7; 0.1, 0.1; 0.5, 0.4 ]; contains three observations from a two-dimensional distribution. While this is the default data arrangement, most functions support different arrangements. It should be noted that the statistics functions don't test for data containing NaN, NA, or Inf. These values need to be detected and dealt with explicitly. See *note isnan: XREFisnan, *note isna: XREFisna, *note isinf: XREFisinf, *note isfinite: XREFisfinite. * Menu: * Descriptive Statistics:: * Basic Statistical Functions:: * Statistical Plots:: * Correlation and Regression Analysis:: * Distributions:: * Tests:: * Random Number Generation::  File: octave.info, Node: Descriptive Statistics, Next: Basic Statistical Functions, Up: Statistics 26.1 Descriptive Statistics =========================== One principal goal of descriptive statistics is to represent the essence of a large data set concisely. Octave provides the mean, median, and mode functions which all summarize a data set with just a single number corresponding to the central tendency of the data. -- Function File: mean (X) -- Function File: mean (X, DIM) -- Function File: mean (X, OPT) -- Function File: mean (X, DIM, OPT) Compute the mean of the elements of the vector X. The mean is defined as mean (x) = SUM_i x(i) / N If X is a matrix, compute the mean for each column and return them in a row vector. If the optional argument DIM is given, operate along this dimension. The optional argument OPT selects the type of mean to compute. The following options are recognized: "a" Compute the (ordinary) arithmetic mean. [default] "g" Compute the geometric mean. "h" Compute the harmonic mean. Both DIM and OPT are optional. If both are supplied, either may appear first. See also: *note median: XREFmedian, *note mode: XREFmode. -- Function File: median (X) -- Function File: median (X, DIM) Compute the median value of the elements of the vector X. When the elements of X are sorted, the median is defined as x(ceil(N/2)) N odd median (x) = (x(N/2) + x((N/2)+1))/2 N even If X is a matrix, compute the median value for each column and return them in a row vector. If the optional DIM argument is given, operate along this dimension. See also: *note mean: XREFmean, *note mode: XREFmode. -- Function File: mode (X) -- Function File: mode (X, DIM) -- Function File: [M, F, C] = mode (...) Compute the most frequently occurring value in a dataset (mode). 'mode' determines the frequency of values along the first non-singleton dimension and returns the value with the highest frequency. If two, or more, values have the same frequency 'mode' returns the smallest. If the optional argument DIM is given, operate along this dimension. The return variable F is the number of occurrences of the mode in the dataset. The cell array C contains all of the elements with the maximum frequency. See also: *note mean: XREFmean, *note median: XREFmedian. Using just one number, such as the mean, to represent an entire data set may not give an accurate picture of the data. One way to characterize the fit is to measure the dispersion of the data. Octave provides several functions for measuring dispersion. -- Function File: range (X) -- Function File: range (X, DIM) Return the range, i.e., the difference between the maximum and the minimum of the input data. If X is a vector, the range is calculated over the elements of X. If X is a matrix, the range is calculated over each column of X. If the optional argument DIM is given, operate along this dimension. The range is a quickly computed measure of the dispersion of a data set, but is less accurate than 'iqr' if there are outlying data points. See also: *note iqr: XREFiqr, *note std: XREFstd. -- Function File: iqr (X) -- Function File: iqr (X, DIM) Return the interquartile range, i.e., the difference between the upper and lower quartile of the input data. If X is a matrix, do the above for first non-singleton dimension of X. If the optional argument DIM is given, operate along this dimension. As a measure of dispersion, the interquartile range is less affected by outliers than either 'range' or 'std'. See also: *note range: XREFrange, *note std: XREFstd. -- Function File: meansq (X) -- Function File: meansq (X, DIM) Compute the mean square of the elements of the vector X. The mean square is defined as meansq (x) = 1/N SUM_i x(i)^2 For matrix arguments, return a row vector containing the mean square of each column. If the optional argument DIM is given, operate along this dimension. See also: *note var: XREFvar, *note std: XREFstd, *note moment: XREFmoment. -- Function File: std (X) -- Function File: std (X, OPT) -- Function File: std (X, OPT, DIM) Compute the standard deviation of the elements of the vector X. The standard deviation is defined as std (x) = sqrt ( 1/(N-1) SUM_i (x(i) - mean(x))^2 ) where N is the number of elements. If X is a matrix, compute the standard deviation for each column and return them in a row vector. The argument OPT determines the type of normalization to use. Valid values are 0: normalize with N-1, provides the square root of the best unbiased estimator of the variance [default] 1: normalize with N, this provides the square root of the second moment around the mean If the optional argument DIM is given, operate along this dimension. See also: *note var: XREFvar, *note range: XREFrange, *note iqr: XREFiqr, *note mean: XREFmean, *note median: XREFmedian. In addition to knowing the size of a dispersion it is useful to know the shape of the data set. For example, are data points massed to the left or right of the mean? Octave provides several common measures to describe the shape of the data set. Octave can also calculate moments allowing arbitrary shape measures to be developed. -- Function File: var (X) -- Function File: var (X, OPT) -- Function File: var (X, OPT, DIM) Compute the variance of the elements of the vector X. The variance is defined as var (x) = 1/(N-1) SUM_i (x(i) - mean(x))^2 If X is a matrix, compute the variance for each column and return them in a row vector. The argument OPT determines the type of normalization to use. Valid values are 0: normalize with N-1, provides the best unbiased estimator of the variance [default] 1: normalizes with N, this provides the second moment around the mean If N==1 the value of OPT is ignored and normalization by N is used. If the optional argument DIM is given, operate along this dimension. See also: *note cov: XREFcov, *note std: XREFstd, *note skewness: XREFskewness, *note kurtosis: XREFkurtosis, *note moment: XREFmoment. -- Function File: skewness (X) -- Function File: skewness (X, FLAG) -- Function File: skewness (X, FLAG, DIM) Compute the sample skewness of the elements of X. The sample skewness is defined as mean ((X - mean (X)).^3) skewness (X) = ------------------------. std (X).^3 The optional argument FLAG controls which normalization is used. If FLAG is equal to 1 (default value, used when FLAG is omitted or empty), return the sample skewness as defined above. If FLAG is equal to 0, return the adjusted skewness coefficient instead: sqrt (N*(N-1)) mean ((X - mean (X)).^3) skewness (X, 0) = -------------- * ------------------------. (N - 2) std (X).^3 The adjusted skewness coefficient is obtained by replacing the sample second and third central moments by their bias-corrected versions. If X is a matrix, or more generally a multi-dimensional array, return the skewness along the first non-singleton dimension. If the optional DIM argument is given, operate along this dimension. See also: *note var: XREFvar, *note kurtosis: XREFkurtosis, *note moment: XREFmoment. -- Function File: kurtosis (X) -- Function File: kurtosis (X, FLAG) -- Function File: kurtosis (X, FLAG, DIM) Compute the sample kurtosis of the elements of X. The sample kurtosis is defined as mean ((X - mean (X)).^4) k1 = ------------------------ std (X).^4 The optional argument FLAG controls which normalization is used. If FLAG is equal to 1 (default value, used when FLAG is omitted or empty), return the sample kurtosis as defined above. If FLAG is equal to 0, return the "bias-corrected" kurtosis coefficient instead: N - 1 k0 = 3 + -------------- * ((N + 1) * k1 - 3 * (N - 1)) (N - 2)(N - 3) The bias-corrected kurtosis coefficient is obtained by replacing the sample second and fourth central moments by their unbiased versions. It is an unbiased estimate of the population kurtosis for normal populations. If X is a matrix, or more generally a multi-dimensional array, return the kurtosis along the first non-singleton dimension. If the optional DIM argument is given, operate along this dimension. See also: *note var: XREFvar, *note skewness: XREFskewness, *note moment: XREFmoment. -- Function File: moment (X, P) -- Function File: moment (X, P, TYPE) -- Function File: moment (X, P, DIM) -- Function File: moment (X, P, TYPE, DIM) -- Function File: moment (X, P, DIM, TYPE) Compute the P-th central moment of the vector X. 1/N SUM_i (x(i) - mean(x))^p If X is a matrix, return the row vector containing the P-th central moment of each column. If the optional argument DIM is given, operate along this dimension. The optional string TYPE specifies the type of moment to be computed. Valid options are: "c" Central Moment (default). "a" "ac" Absolute Central Moment. The moment about the mean ignoring sign defined as 1/N SUM_i (abs (x(i) - mean(x)))^p "r" Raw Moment. The moment about zero defined as moment (x) = 1/N SUM_i x(i)^p "ar" Absolute Raw Moment. The moment about zero ignoring sign defined as 1/N SUM_i ( abs (x(i)) )^p If both TYPE and DIM are given they may appear in any order. See also: *note var: XREFvar, *note skewness: XREFskewness, *note kurtosis: XREFkurtosis. -- Function File: Q = quantile (X) -- Function File: Q = quantile (X, P) -- Function File: Q = quantile (X, P, DIM) -- Function File: Q = quantile (X, P, DIM, METHOD) For a sample, X, calculate the quantiles, Q, corresponding to the cumulative probability values in P. All non-numeric values (NaNs) of X are ignored. If X is a matrix, compute the quantiles for each column and return them in a matrix, such that the i-th row of Q contains the P(i)th quantiles of each column of X. If P is unspecified, return the quantiles for '[0.00 0.25 0.50 0.75 1.00]'. The optional argument DIM determines the dimension along which the quantiles are calculated. If DIM is omitted it defaults to the first non-singleton dimension. The methods available to calculate sample quantiles are the nine methods used by R (). The default value is METHOD = 5. Discontinuous sample quantile methods 1, 2, and 3 1. Method 1: Inverse of empirical distribution function. 2. Method 2: Similar to method 1 but with averaging at discontinuities. 3. Method 3: SAS definition: nearest even order statistic. Continuous sample quantile methods 4 through 9, where p(k) is the linear interpolation function respecting each methods' representative cdf. 4. Method 4: p(k) = k / n. That is, linear interpolation of the empirical cdf. 5. Method 5: p(k) = (k - 0.5) / n. That is a piecewise linear function where the knots are the values midway through the steps of the empirical cdf. 6. Method 6: p(k) = k / (n + 1). 7. Method 7: p(k) = (k - 1) / (n - 1). 8. Method 8: p(k) = (k - 1/3) / (n + 1/3). The resulting quantile estimates are approximately median-unbiased regardless of the distribution of X. 9. Method 9: p(k) = (k - 3/8) / (n + 1/4). The resulting quantile estimates are approximately unbiased for the expected order statistics if X is normally distributed. Hyndman and Fan (1996) recommend method 8. Maxima, S, and R (versions prior to 2.0.0) use 7 as their default. Minitab and SPSS use method 6. MATLAB uses method 5. References: * Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole. * Hyndman, R. J. and Fan, Y. (1996) Sample quantiles in statistical packages, American Statistician, 50, 361-365. * R: A Language and Environment for Statistical Computing; . Examples: x = randi (1000, [10, 1]); # Create empirical data in range 1-1000 q = quantile (x, [0, 1]); # Return minimum, maximum of distribution q = quantile (x, [0.25 0.5 0.75]); # Return quartiles of distribution See also: *note prctile: XREFprctile. -- Function File: Q = prctile (X) -- Function File: Q = prctile (X, P) -- Function File: Q = prctile (X, P, DIM) For a sample X, compute the quantiles, Q, corresponding to the cumulative probability values, P, in percent. If X is a matrix, compute the percentiles for each column and return them in a matrix, such that the i-th row of Y contains the P(i)th percentiles of each column of X. If P is unspecified, return the quantiles for '[0 25 50 75 100]'. The optional argument DIM determines the dimension along which the percentiles are calculated. If DIM is omitted it defaults to the first non-singleton dimension. Programming Note: All non-numeric values (NaNs) of X are ignored. See also: *note quantile: XREFquantile. A summary view of a data set can be generated quickly with the 'statistics' function. -- Function File: statistics (X) -- Function File: statistics (X, DIM) Return a vector with the minimum, first quartile, median, third quartile, maximum, mean, standard deviation, skewness, and kurtosis of the elements of the vector X. If X is a matrix, calculate statistics over the first non-singleton dimension. If the optional argument DIM is given, operate along this dimension. See also: *note min: XREFmin, *note max: XREFmax, *note median: XREFmedian, *note mean: XREFmean, *note std: XREFstd, *note skewness: XREFskewness, *note kurtosis: XREFkurtosis.  File: octave.info, Node: Basic Statistical Functions, Next: Statistical Plots, Prev: Descriptive Statistics, Up: Statistics 26.2 Basic Statistical Functions ================================ Octave supports various helpful statistical functions. Many are useful as initial steps to prepare a data set for further analysis. Others provide different measures from those of the basic descriptive statistics. -- Function File: center (X) -- Function File: center (X, DIM) Center data by subtracting its mean. If X is a vector, subtract its mean. If X is a matrix, do the above for each column. If the optional argument DIM is given, operate along this dimension. Programming Note: 'center' has obvious application for normalizing statistical data. It is also useful for improving the precision of general numerical calculations. Whenever there is a large value that is common to a batch of data, the mean can be subtracted off, the calculation performed, and then the mean added back to obtain the final answer. See also: *note zscore: XREFzscore. -- Function File: Z = zscore (X) -- Function File: Z = zscore (X, OPT) -- Function File: Z = zscore (X, OPT, DIM) -- Function File: [Z, MU, SIGMA] = zscore (...) Compute the Z score of X If X is a vector, subtract its mean and divide by its standard deviation. If the standard deviation is zero, divide by 1 instead. The optional parameter OPT determines the normalization to use when computing the standard deviation and has the same definition as the corresponding parameter for 'std'. If X is a matrix, calculate along the first non-singleton dimension. If the third optional argument DIM is given, operate along this dimension. The optional outputs MU and SIGMA contain the mean and standard deviation. See also: *note mean: XREFmean, *note std: XREFstd, *note center: XREFcenter. -- Function File: N = histc (X, EDGES) -- Function File: N = histc (X, EDGES, DIM) -- Function File: [N, IDX] = histc (...) Compute histogram counts. When X is a vector, the function counts the number of elements of X that fall in the histogram bins defined by EDGES. This must be a vector of monotonically increasing values that define the edges of the histogram bins. 'N(k)' contains the number of elements in X for which 'EDGES(k) <= X < EDGES(k+1)'. The final element of N contains the number of elements of X exactly equal to the last element of EDGES. When X is an N-dimensional array, the computation is carried out along dimension DIM. If not specified DIM defaults to the first non-singleton dimension. When a second output argument is requested an index matrix is also returned. The IDX matrix has the same size as X. Each element of IDX contains the index of the histogram bin in which the corresponding element of X was counted. See also: *note hist: XREFhist. -- Function File: C = nchoosek (N, K) -- Function File: C = nchoosek (SET, K) Compute the binomial coefficient of N or list all possible combinations of a SET of items. If N is a scalar then calculate the binomial coefficient of N and K which is defined as / \ | n | n (n-1) (n-2) ... (n-k+1) n! | | = ------------------------- = --------- | k | k! k! (n-k)! \ / This is the number of combinations of N items taken in groups of size K. If the first argument is a vector, SET, then generate all combinations of the elements of SET, taken K at a time, with one row per combination. The result C has K columns and 'nchoosek (length (SET), K)' rows. For example: How many ways can three items be grouped into pairs? nchoosek (3, 2) => 3 What are the possible pairs? nchoosek (1:3, 2) => 1 2 1 3 2 3 Programming Note: When calculating the binomial coefficient 'nchoosek' works only for non-negative, integer arguments. Use 'bincoeff' for non-integer and negative scalar arguments, or for computing many binomial coefficients at once with vector inputs for N or K. See also: *note bincoeff: XREFbincoeff, *note perms: XREFperms. -- Function File: perms (V) Generate all permutations of V with one row per permutation. The result has size 'factorial (N) * N', where N is the length of V. Example perms ([1, 2, 3]) => 1 2 3 2 1 3 1 3 2 2 3 1 3 1 2 3 2 1 Programming Note: The maximum length of V should be less than or equal to 10 to limit memory consumption. See also: *note permute: XREFpermute, *note randperm: XREFrandperm, *note nchoosek: XREFnchoosek. -- Function File: ranks (X, DIM) Return the ranks of X along the first non-singleton dimension adjusted for ties. If the optional argument DIM is given, operate along this dimension. See also: *note spearman: XREFspearman, *note kendall: XREFkendall. -- Function File: run_count (X, N) -- Function File: run_count (X, N, DIM) Count the upward runs along the first non-singleton dimension of X of length 1, 2, ..., N-1 and greater than or equal to N. If the optional argument DIM is given then operate along this dimension. See also: *note runlength: XREFrunlength. -- Function File: count = runlength (X) -- Function File: [count, value] = runlength (X) Find the lengths of all sequences of common values. COUNT is a vector with the lengths of each repeated value. The optional output VALUE contains the value that was repeated in the sequence. runlength ([2, 2, 0, 4, 4, 4, 0, 1, 1, 1, 1]) => [2, 1, 3, 1, 4] See also: *note run_count: XREFrun_count. -- Function File: probit (P) Return the probit (the quantile of the standard normal distribution) for each element of P. See also: *note logit: XREFlogit. -- Function File: logit (P) Compute the logit for each value of P The logit is defined as logit (P) = log (P / (1-P)) See also: *note probit: XREFprobit, *note logistic_cdf: XREFlogistic_cdf. -- Function File: cloglog (X) Return the complementary log-log function of X. The complementary log-log function is defined as cloglog (x) = - log (- log (X)) -- Function File: mahalanobis (X, Y) Return the Mahalanobis' D-square distance between the multivariate samples X and Y. The data X and Y must have the same number of components (columns), but may have a different number of observations (rows). -- Function File: [T, L_X] = table (X) -- Function File: [T, L_X, L_Y] = table (X, Y) Create a contingency table T from data vectors. The L_X and L_Y vectors are the corresponding levels. Currently, only 1- and 2-dimensional tables are supported.  File: octave.info, Node: Statistical Plots, Next: Correlation and Regression Analysis, Prev: Basic Statistical Functions, Up: Statistics 26.3 Statistical Plots ====================== Octave can create Quantile Plots (QQ-Plots), and Probability Plots (PP-Plots). These are simple graphical tests for determining if a data set comes from a certain distribution. Note that Octave can also show histograms of data using the 'hist' function as described in *note Two-Dimensional Plots::. -- Function File: [Q, S] = qqplot (X) -- Function File: [Q, S] = qqplot (X, Y) -- Function File: [Q, S] = qqplot (X, DIST) -- Function File: [Q, S] = qqplot (X, Y, PARAMS) -- Function File: qqplot (...) Perform a QQ-plot (quantile plot). If F is the CDF of the distribution DIST with parameters PARAMS and G its inverse, and X a sample vector of length N, the QQ-plot graphs ordinate S(I) = I-th largest element of x versus abscissa Q(If) = G((I - 0.5)/N). If the sample comes from F, except for a transformation of location and scale, the pairs will approximately follow a straight line. If the second argument is a vector Y the empirical CDF of Y is used as DIST. The default for DIST is the standard normal distribution. The optional argument PARAMS contains a list of parameters of DIST. For example, for a quantile plot of the uniform distribution on [2,4] and X, use qqplot (x, "unif", 2, 4) DIST can be any string for which a function DISTINV or DIST_INV exists that calculates the inverse CDF of distribution DIST. If no output arguments are given, the data are plotted directly. -- Function File: [P, Y] = ppplot (X, DIST, PARAMS) Perform a PP-plot (probability plot). If F is the CDF of the distribution DIST with parameters PARAMS and X a sample vector of length N, the PP-plot graphs ordinate Y(I) = F (I-th largest element of X) versus abscissa P(I) = (I - 0.5)/N. If the sample comes from F, the pairs will approximately follow a straight line. The default for DIST is the standard normal distribution. The optional argument PARAMS contains a list of parameters of DIST. For example, for a probability plot of the uniform distribution on [2,4] and X, use ppplot (x, "uniform", 2, 4) DIST can be any string for which a function DIST_CDF that calculates the CDF of distribution DIST exists. If no output is requested then the data are plotted immediately.  File: octave.info, Node: Correlation and Regression Analysis, Next: Distributions, Prev: Statistical Plots, Up: Statistics 26.4 Correlation and Regression Analysis ======================================== -- Function File: cov (X) -- Function File: cov (X, OPT) -- Function File: cov (X, Y) -- Function File: cov (X, Y, OPT) Compute the covariance matrix. If each row of X and Y is an observation, and each column is a variable, then the (I, J)-th entry of 'cov (X, Y)' is the covariance between the I-th variable in X and the J-th variable in Y. cov (x) = 1/N-1 * SUM_i (x(i) - mean(x)) * (y(i) - mean(y)) If called with one argument, compute 'cov (X, X)', the covariance between the columns of X. The argument OPT determines the type of normalization to use. Valid values are 0: normalize with N-1, provides the best unbiased estimator of the covariance [default] 1: normalize with N, this provides the second moment around the mean Compatibility Note:: Octave always computes the covariance matrix. For two inputs, however, MATLAB will calculate 'cov (X(:), Y(:))' whenever the number of elements in X and Y are equal. This will result in a scalar rather than a matrix output. Code relying on this odd definition will need to be changed when running in Octave. See also: *note corr: XREFcorr. -- Function File: corr (X) -- Function File: corr (X, Y) Compute matrix of correlation coefficients. If each row of X and Y is an observation and each column is a variable, then the (I, J)-th entry of 'corr (X, Y)' is the correlation between the I-th variable in X and the J-th variable in Y. corr (x,y) = cov (x,y) / (std (x) * std (y)) If called with one argument, compute 'corr (X, X)', the correlation between the columns of X. See also: *note cov: XREFcov. -- Function File: spearman (X) -- Function File: spearman (X, Y) Compute Spearman's rank correlation coefficient RHO. For two data vectors X and Y, Spearman's RHO is the correlation coefficient of the ranks of X and Y. If X and Y are drawn from independent distributions, RHO has zero mean and variance '1 / (n - 1)', and is asymptotically normally distributed. 'spearman (X)' is equivalent to 'spearman (X, X)'. See also: *note ranks: XREFranks, *note kendall: XREFkendall. -- Function File: kendall (X) -- Function File: kendall (X, Y) Compute Kendall's TAU. For two data vectors X, Y of common length N, Kendall's TAU is the correlation of the signs of all rank differences of X and Y; i.e., if both X and Y have distinct entries, then 1 tau = ------- SUM sign (q(i) - q(j)) * sign (r(i) - r(j)) n (n-1) i,j in which the Q(I) and R(I) are the ranks of X and Y, respectively. If X and Y are drawn from independent distributions, Kendall's TAU is asymptotically normal with mean 0 and variance '(2 * (2N+5)) / (9 * N * (N-1))'. 'kendall (X)' is equivalent to 'kendall (X, X)'. See also: *note ranks: XREFranks, *note spearman: XREFspearman. -- Function File: [THETA, BETA, DEV, DL, D2L, P] = logistic_regression (Y, X, PRINT, THETA, BETA) Perform ordinal logistic regression. Suppose Y takes values in K ordered categories, and let 'gamma_i (X)' be the cumulative probability that Y falls in one of the first I categories given the covariate X. Then [theta, beta] = logistic_regression (y, x) fits the model logit (gamma_i (x)) = theta_i - beta' * x, i = 1 ... k-1 The number of ordinal categories, K, is taken to be the number of distinct values of 'round (Y)'. If K equals 2, Y is binary and the model is ordinary logistic regression. The matrix X is assumed to have full column rank. Given Y only, 'theta = logistic_regression (y)' fits the model with baseline logit odds only. The full form is [theta, beta, dev, dl, d2l, gamma] = logistic_regression (y, x, print, theta, beta) in which all output arguments and all input arguments except Y are optional. Setting PRINT to 1 requests summary information about the fitted model to be displayed. Setting PRINT to 2 requests information about convergence at each iteration. Other values request no information to be displayed. The input arguments THETA and BETA give initial estimates for THETA and BETA. The returned value DEV holds minus twice the log-likelihood. The returned values DL and D2L are the vector of first and the matrix of second derivatives of the log-likelihood with respect to THETA and BETA. P holds estimates for the conditional distribution of Y given X.  File: octave.info, Node: Distributions, Next: Tests, Prev: Correlation and Regression Analysis, Up: Statistics 26.5 Distributions ================== Octave has functions for computing the Probability Density Function (PDF), the Cumulative Distribution function (CDF), and the quantile (the inverse of the CDF) for a large number of distributions. The following table summarizes the supported distributions (in alphabetical order). Distribution PDF CDF Quantile ----------------------------------------------------------------------------- Beta Distribution 'betapdf' 'betacdf' 'betainv' Binomial 'binopdf' 'binocdf' 'binoinv' Distribution Cauchy Distribution 'cauchy_pdf' 'cauchy_cdf' 'cauchy_inv' Chi-Square 'chi2pdf' 'chi2cdf' 'chi2inv' Distribution Univariate Discrete 'discrete_pdf' 'discrete_cdf' 'discrete_inv' Distribution Empirical 'empirical_pdf' 'empirical_cdf' 'empirical_inv' Distribution Exponential 'exppdf' 'expcdf' 'expinv' Distribution F Distribution 'fpdf' 'fcdf' 'finv' Gamma Distribution 'gampdf' 'gamcdf' 'gaminv' Geometric 'geopdf' 'geocdf' 'geoinv' Distribution Hypergeometric 'hygepdf' 'hygecdf' 'hygeinv' Distribution Kolmogorov Smirnov _Not Available_ 'kolmogorov_smirnov_cdf'_Not Available_ Distribution Laplace Distribution 'laplace_pdf' 'laplace_cdf' 'laplace_inv' Logistic 'logistic_pdf' 'logistic_cdf' 'logistic_inv' Distribution Log-Normal 'lognpdf' 'logncdf' 'logninv' Distribution Univariate Normal 'normpdf' 'normcdf' 'norminv' Distribution Pascal Distribution 'nbinpdf' 'nbincdf' 'nbininv' Poisson Distribution 'poisspdf' 'poisscdf' 'poissinv' Standard Normal 'stdnormal_pdf' 'stdnormal_cdf' 'stdnormal_inv' Distribution t (Student) 'tpdf' 'tcdf' 'tinv' Distribution Univariate Discrete 'unidpdf' 'unidcdf' 'unidinv' Distribution Uniform Distribution 'unifpdf' 'unifcdf' 'unifinv' Weibull Distribution 'wblpdf' 'wblcdf' 'wblinv' -- Function File: betapdf (X, A, B) For each element of X, compute the probability density function (PDF) at X of the Beta distribution with parameters A and B. -- Function File: betacdf (X, A, B) For each element of X, compute the cumulative distribution function (CDF) at X of the Beta distribution with parameters A and B. -- Function File: betainv (X, A, B) For each element of X, compute the quantile (the inverse of the CDF) at X of the Beta distribution with parameters A and B. -- Function File: binopdf (X, N, P) For each element of X, compute the probability density function (PDF) at X of the binomial distribution with parameters N and P, where N is the number of trials and P is the probability of success. -- Function File: binocdf (X, N, P) For each element of X, compute the cumulative distribution function (CDF) at X of the binomial distribution with parameters N and P, where N is the number of trials and P is the probability of success. -- Function File: binoinv (X, N, P) For each element of X, compute the quantile (the inverse of the CDF) at X of the binomial distribution with parameters N and P, where N is the number of trials and P is the probability of success. -- Function File: cauchy_pdf (X) -- Function File: cauchy_pdf (X, LOCATION, SCALE) For each element of X, compute the probability density function (PDF) at X of the Cauchy distribution with location parameter LOCATION and scale parameter SCALE > 0. Default values are LOCATION = 0, SCALE = 1. -- Function File: cauchy_cdf (X) -- Function File: cauchy_cdf (X, LOCATION, SCALE) For each element of X, compute the cumulative distribution function (CDF) at X of the Cauchy distribution with location parameter LOCATION and scale parameter SCALE. Default values are LOCATION = 0, SCALE = 1. -- Function File: cauchy_inv (X) -- Function File: cauchy_inv (X, LOCATION, SCALE) For each element of X, compute the quantile (the inverse of the CDF) at X of the Cauchy distribution with location parameter LOCATION and scale parameter SCALE. Default values are LOCATION = 0, SCALE = 1. -- Function File: chi2pdf (X, N) For each element of X, compute the probability density function (PDF) at X of the chi-square distribution with N degrees of freedom. -- Function File: chi2cdf (X, N) For each element of X, compute the cumulative distribution function (CDF) at X of the chi-square distribution with N degrees of freedom. -- Function File: chi2inv (X, N) For each element of X, compute the quantile (the inverse of the CDF) at X of the chi-square distribution with N degrees of freedom. -- Function File: discrete_pdf (X, V, P) For each element of X, compute the probability density function (PDF) at X of a univariate discrete distribution which assumes the values in V with probabilities P. -- Function File: discrete_cdf (X, V, P) For each element of X, compute the cumulative distribution function (CDF) at X of a univariate discrete distribution which assumes the values in V with probabilities P. -- Function File: discrete_inv (X, V, P) For each element of X, compute the quantile (the inverse of the CDF) at X of the univariate distribution which assumes the values in V with probabilities P. -- Function File: empirical_pdf (X, DATA) For each element of X, compute the probability density function (PDF) at X of the empirical distribution obtained from the univariate sample DATA. -- Function File: empirical_cdf (X, DATA) For each element of X, compute the cumulative distribution function (CDF) at X of the empirical distribution obtained from the univariate sample DATA. -- Function File: empirical_inv (X, DATA) For each element of X, compute the quantile (the inverse of the CDF) at X of the empirical distribution obtained from the univariate sample DATA. -- Function File: exppdf (X, LAMBDA) For each element of X, compute the probability density function (PDF) at X of the exponential distribution with mean LAMBDA. -- Function File: expcdf (X, LAMBDA) For each element of X, compute the cumulative distribution function (CDF) at X of the exponential distribution with mean LAMBDA. The arguments can be of common size or scalars. -- Function File: expinv (X, LAMBDA) For each element of X, compute the quantile (the inverse of the CDF) at X of the exponential distribution with mean LAMBDA. -- Function File: fpdf (X, M, N) For each element of X, compute the probability density function (PDF) at X of the F distribution with M and N degrees of freedom. -- Function File: fcdf (X, M, N) For each element of X, compute the cumulative distribution function (CDF) at X of the F distribution with M and N degrees of freedom. -- Function File: finv (X, M, N) For each element of X, compute the quantile (the inverse of the CDF) at X of the F distribution with M and N degrees of freedom. -- Function File: gampdf (X, A, B) For each element of X, return the probability density function (PDF) at X of the Gamma distribution with shape parameter A and scale B. -- Function File: gamcdf (X, A, B) For each element of X, compute the cumulative distribution function (CDF) at X of the Gamma distribution with shape parameter A and scale B. -- Function File: gaminv (X, A, B) For each element of X, compute the quantile (the inverse of the CDF) at X of the Gamma distribution with shape parameter A and scale B. -- Function File: geopdf (X, P) For each element of X, compute the probability density function (PDF) at X of the geometric distribution with parameter P. The geometric distribution models the number of failures (X-1) of a Bernoulli trial with probability P before the first success (X). -- Function File: geocdf (X, P) For each element of X, compute the cumulative distribution function (CDF) at X of the geometric distribution with parameter P. The geometric distribution models the number of failures (X-1) of a Bernoulli trial with probability P before the first success (X). -- Function File: geoinv (X, P) For each element of X, compute the quantile (the inverse of the CDF) at X of the geometric distribution with parameter P. The geometric distribution models the number of failures (X-1) of a Bernoulli trial with probability P before the first success (X). -- Function File: hygepdf (X, T, M, N) Compute the probability density function (PDF) at X of the hypergeometric distribution with parameters T, M, and N. This is the probability of obtaining X marked items when randomly drawing a sample of size N without replacement from a population of total size T containing M marked items. The parameters T, M, and N must be positive integers with M and N not greater than T. -- Function File: hygecdf (X, T, M, N) Compute the cumulative distribution function (CDF) at X of the hypergeometric distribution with parameters T, M, and N. This is the probability of obtaining not more than X marked items when randomly drawing a sample of size N without replacement from a population of total size T containing M marked items. The parameters T, M, and N must be positive integers with M and N not greater than T. -- Function File: hygeinv (X, T, M, N) For each element of X, compute the quantile (the inverse of the CDF) at X of the hypergeometric distribution with parameters T, M, and N. This is the probability of obtaining X marked items when randomly drawing a sample of size N without replacement from a population of total size T containing M marked items. The parameters T, M, and N must be positive integers with M and N not greater than T. -- Function File: kolmogorov_smirnov_cdf (X, TOL) Return the cumulative distribution function (CDF) at X of the Kolmogorov-Smirnov distribution. This is defined as Inf Q(x) = SUM (-1)^k exp (-2 k^2 x^2) k = -Inf for X > 0. The optional parameter TOL specifies the precision up to which the series should be evaluated; the default is TOL = 'eps'. -- Function File: laplace_pdf (X) For each element of X, compute the probability density function (PDF) at X of the Laplace distribution. -- Function File: laplace_cdf (X) For each element of X, compute the cumulative distribution function (CDF) at X of the Laplace distribution. -- Function File: laplace_inv (X) For each element of X, compute the quantile (the inverse of the CDF) at X of the Laplace distribution. -- Function File: logistic_pdf (X) For each element of X, compute the PDF at X of the logistic distribution. -- Function File: logistic_cdf (X) For each element of X, compute the cumulative distribution function (CDF) at X of the logistic distribution. -- Function File: logistic_inv (X) For each element of X, compute the quantile (the inverse of the CDF) at X of the logistic distribution. -- Function File: lognpdf (X) -- Function File: lognpdf (X, MU, SIGMA) For each element of X, compute the probability density function (PDF) at X of the lognormal distribution with parameters MU and SIGMA. If a random variable follows this distribution, its logarithm is normally distributed with mean MU and standard deviation SIGMA. Default values are MU = 0, SIGMA = 1. -- Function File: logncdf (X) -- Function File: logncdf (X, MU, SIGMA) For each element of X, compute the cumulative distribution function (CDF) at X of the lognormal distribution with parameters MU and SIGMA. If a random variable follows this distribution, its logarithm is normally distributed with mean MU and standard deviation SIGMA. Default values are MU = 0, SIGMA = 1. -- Function File: logninv (X) -- Function File: logninv (X, MU, SIGMA) For each element of X, compute the quantile (the inverse of the CDF) at X of the lognormal distribution with parameters MU and SIGMA. If a random variable follows this distribution, its logarithm is normally distributed with mean MU and standard deviation SIGMA. Default values are MU = 0, SIGMA = 1. -- Function File: nbinpdf (X, N, P) For each element of X, compute the probability density function (PDF) at X of the negative binomial distribution with parameters N and P. When N is integer this is the Pascal distribution. When N is extended to real numbers this is the Polya distribution. The number of failures in a Bernoulli experiment with success probability P before the N-th success follows this distribution. -- Function File: nbincdf (X, N, P) For each element of X, compute the cumulative distribution function (CDF) at X of the negative binomial distribution with parameters N and P. When N is integer this is the Pascal distribution. When N is extended to real numbers this is the Polya distribution. The number of failures in a Bernoulli experiment with success probability P before the N-th success follows this distribution. -- Function File: nbininv (X, N, P) For each element of X, compute the quantile (the inverse of the CDF) at X of the negative binomial distribution with parameters N and P. When N is integer this is the Pascal distribution. When N is extended to real numbers this is the Polya distribution. The number of failures in a Bernoulli experiment with success probability P before the N-th success follows this distribution. -- Function File: normpdf (X) -- Function File: normpdf (X, MU, SIGMA) For each element of X, compute the probability density function (PDF) at X of the normal distribution with mean MU and standard deviation SIGMA. Default values are MU = 0, SIGMA = 1. -- Function File: normcdf (X) -- Function File: normcdf (X, MU, SIGMA) For each element of X, compute the cumulative distribution function (CDF) at X of the normal distribution with mean MU and standard deviation SIGMA. Default values are MU = 0, SIGMA = 1. -- Function File: norminv (X) -- Function File: norminv (X, MU, SIGMA) For each element of X, compute the quantile (the inverse of the CDF) at X of the normal distribution with mean MU and standard deviation SIGMA. Default values are MU = 0, SIGMA = 1. -- Function File: poisspdf (X, LAMBDA) For each element of X, compute the probability density function (PDF) at X of the Poisson distribution with parameter LAMBDA. -- Function File: poisscdf (X, LAMBDA) For each element of X, compute the cumulative distribution function (CDF) at X of the Poisson distribution with parameter LAMBDA. -- Function File: poissinv (X, LAMBDA) For each element of X, compute the quantile (the inverse of the CDF) at X of the Poisson distribution with parameter LAMBDA. -- Function File: stdnormal_pdf (X) For each element of X, compute the probability density function (PDF) at X of the standard normal distribution (mean = 0, standard deviation = 1). -- Function File: stdnormal_cdf (X) For each element of X, compute the cumulative distribution function (CDF) at X of the standard normal distribution (mean = 0, standard deviation = 1). -- Function File: stdnormal_inv (X) For each element of X, compute the quantile (the inverse of the CDF) at X of the standard normal distribution (mean = 0, standard deviation = 1). -- Function File: tpdf (X, N) For each element of X, compute the probability density function (PDF) at X of the T (Student) distribution with N degrees of freedom. -- Function File: tcdf (X, N) For each element of X, compute the cumulative distribution function (CDF) at X of the t (Student) distribution with N degrees of freedom. -- Function File: tinv (X, N) For each element of X, compute the quantile (the inverse of the CDF) at X of the t (Student) distribution with N degrees of freedom. This function is analogous to looking in a table for the t-value of a single-tailed distribution. -- Function File: unidpdf (X, N) For each element of X, compute the probability density function (PDF) at X of a discrete uniform distribution which assumes the integer values 1-N with equal probability. Warning: The underlying implementation uses the double class and will only be accurate for N <= 'bitmax' (2^{53} - 1 on IEEE-754 compatible systems). -- Function File: unidcdf (X, N) For each element of X, compute the cumulative distribution function (CDF) at X of a discrete uniform distribution which assumes the integer values 1-N with equal probability. -- Function File: unidinv (X, N) For each element of X, compute the quantile (the inverse of the CDF) at X of the discrete uniform distribution which assumes the integer values 1-N with equal probability. -- Function File: unifpdf (X) -- Function File: unifpdf (X, A, B) For each element of X, compute the probability density function (PDF) at X of the uniform distribution on the interval [A, B]. Default values are A = 0, B = 1. -- Function File: unifcdf (X) -- Function File: unifcdf (X, A, B) For each element of X, compute the cumulative distribution function (CDF) at X of the uniform distribution on the interval [A, B]. Default values are A = 0, B = 1. -- Function File: unifinv (X) -- Function File: unifinv (X, A, B) For each element of X, compute the quantile (the inverse of the CDF) at X of the uniform distribution on the interval [A, B]. Default values are A = 0, B = 1. -- Function File: wblpdf (X) -- Function File: wblpdf (X, SCALE) -- Function File: wblpdf (X, SCALE, SHAPE) Compute the probability density function (PDF) at X of the Weibull distribution with scale parameter SCALE and shape parameter SHAPE. This is given by shape * scale^(-shape) * x^(shape-1) * exp (-(x/scale)^shape) for X >= 0. Default values are SCALE = 1, SHAPE = 1. -- Function File: wblcdf (X) -- Function File: wblcdf (X, SCALE) -- Function File: wblcdf (X, SCALE, SHAPE) Compute the cumulative distribution function (CDF) at X of the Weibull distribution with scale parameter SCALE and shape parameter SHAPE. This is defined as 1 - exp (-(x/scale)^shape) for X >= 0. Default values are SCALE = 1, SHAPE = 1. -- Function File: wblinv (X) -- Function File: wblinv (X, SCALE) -- Function File: wblinv (X, SCALE, SHAPE) Compute the quantile (the inverse of the CDF) at X of the Weibull distribution with scale parameter SCALE and shape parameter SHAPE. Default values are SCALE = 1, SHAPE = 1.  File: octave.info, Node: Tests, Next: Random Number Generation, Prev: Distributions, Up: Statistics 26.6 Tests ========== Octave can perform many different statistical tests. The following table summarizes the available tests. Hypothesis Test Functions ------------------------------------------------------------------- Equal mean values 'anova', 'hotelling_test2', 't_test_2', 'welch_test', 'wilcoxon_test', 'z_test_2' Equal medians 'kruskal_wallis_test', 'sign_test' Equal variances 'bartlett_test', 'manova', 'var_test' Equal distributions 'chisquare_test_homogeneity', 'kolmogorov_smirnov_test_2', 'u_test' Equal marginal frequencies 'mcnemar_test' Equal success probabilities 'prop_test_2' Independent observations 'chisquare_test_independence', 'run_test' Uncorrelated observations 'cor_test' Given mean value 'hotelling_test', 't_test', 'z_test' Observations from given 'kolmogorov_smirnov_test' distribution Regression 'f_test_regression', 't_test_regression' The tests return a p-value that describes the outcome of the test. Assuming that the test hypothesis is true, the p-value is the probability of obtaining a worse result than the observed one. So large p-values corresponds to a successful test. Usually a test hypothesis is accepted if the p-value exceeds 0.05. -- Function File: [PVAL, F, DF_B, DF_W] = anova (Y, G) Perform a one-way analysis of variance (ANOVA). The goal is to test whether the population means of data taken from K different groups are all equal. Data may be given in a single vector Y with groups specified by a corresponding vector of group labels G (e.g., numbers from 1 to K). This is the general form which does not impose any restriction on the number of data in each group or the group labels. If Y is a matrix and G is omitted, each column of Y is treated as a group. This form is only appropriate for balanced ANOVA in which the numbers of samples from each group are all equal. Under the null of constant means, the statistic F follows an F distribution with DF_B and DF_W degrees of freedom. The p-value (1 minus the CDF of this distribution at F) is returned in PVAL. If no output argument is given, the standard one-way ANOVA table is printed. See also: *note manova: XREFmanova. -- Function File: [PVAL, CHISQ, DF] = bartlett_test (X1, ...) Perform a Bartlett test for the homogeneity of variances in the data vectors X1, X2, ..., XK, where K > 1. Under the null of equal variances, the test statistic CHISQ approximately follows a chi-square distribution with DF degrees of freedom. The p-value (1 minus the CDF of this distribution at CHISQ) is returned in PVAL. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, CHISQ, DF] = chisquare_test_homogeneity (X, Y, C) Given two samples X and Y, perform a chisquare test for homogeneity of the null hypothesis that X and Y come from the same distribution, based on the partition induced by the (strictly increasing) entries of C. For large samples, the test statistic CHISQ approximately follows a chisquare distribution with DF = 'length (C)' degrees of freedom. The p-value (1 minus the CDF of this distribution at CHISQ) is returned in PVAL. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, CHISQ, DF] = chisquare_test_independence (X) Perform a chi-square test for independence based on the contingency table X. Under the null hypothesis of independence, CHISQ approximately has a chi-square distribution with DF degrees of freedom. The p-value (1 minus the CDF of this distribution at chisq) of the test is returned in PVAL. If no output argument is given, the p-value is displayed. -- Function File: cor_test (X, Y, ALT, METHOD) Test whether two samples X and Y come from uncorrelated populations. The optional argument string ALT describes the alternative hypothesis, and can be "!=" or "<>" (nonzero), ">" (greater than 0), or "<" (less than 0). The default is the two-sided case. The optional argument string METHOD specifies which correlation coefficient to use for testing. If METHOD is "pearson" (default), the (usual) Pearson's produt moment correlation coefficient is used. In this case, the data should come from a bivariate normal distribution. Otherwise, the other two methods offer nonparametric alternatives. If METHOD is "kendall", then Kendall's rank correlation tau is used. If METHOD is "spearman", then Spearman's rank correlation rho is used. Only the first character is necessary. The output is a structure with the following elements: PVAL The p-value of the test. STAT The value of the test statistic. DIST The distribution of the test statistic. PARAMS The parameters of the null distribution of the test statistic. ALTERNATIVE The alternative hypothesis. METHOD The method used for testing. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, F, DF_NUM, DF_DEN] = f_test_regression (Y, X, RR, R) Perform an F test for the null hypothesis rr * b = r in a classical normal regression model y = X * b + e. Under the null, the test statistic F follows an F distribution with DF_NUM and DF_DEN degrees of freedom. The p-value (1 minus the CDF of this distribution at F) is returned in PVAL. If not given explicitly, R = 0. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, TSQ] = hotelling_test (X, M) For a sample X from a multivariate normal distribution with unknown mean and covariance matrix, test the null hypothesis that 'mean (X) == M'. Hotelling's T^2 is returned in TSQ. Under the null, (n-p) T^2 / (p(n-1)) has an F distribution with p and n-p degrees of freedom, where n and p are the numbers of samples and variables, respectively. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, TSQ] = hotelling_test_2 (X, Y) For two samples X from multivariate normal distributions with the same number of variables (columns), unknown means and unknown equal covariance matrices, test the null hypothesis 'mean (X) == mean (Y)'. Hotelling's two-sample T^2 is returned in TSQ. Under the null, (n_x+n_y-p-1) T^2 / (p(n_x+n_y-2)) has an F distribution with p and n_x+n_y-p-1 degrees of freedom, where n_x and n_y are the sample sizes and p is the number of variables. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, KS] = kolmogorov_smirnov_test (X, DIST, PARAMS, ALT) Perform a Kolmogorov-Smirnov test of the null hypothesis that the sample X comes from the (continuous) distribution DIST. if F and G are the CDFs corresponding to the sample and dist, respectively, then the null is that F == G. The optional argument PARAMS contains a list of parameters of DIST. For example, to test whether a sample X comes from a uniform distribution on [2,4], use kolmogorov_smirnov_test (x, "unif", 2, 4) DIST can be any string for which a function DISTCDF that calculates the CDF of distribution DIST exists. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative F != G. In this case, the test statistic KS follows a two-sided Kolmogorov-Smirnov distribution. If ALT is ">", the one-sided alternative F > G is considered. Similarly for "<", the one-sided alternative F > G is considered. In this case, the test statistic KS has a one-sided Kolmogorov-Smirnov distribution. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, KS, D] = kolmogorov_smirnov_test_2 (X, Y, ALT) Perform a 2-sample Kolmogorov-Smirnov test of the null hypothesis that the samples X and Y come from the same (continuous) distribution. If F and G are the CDFs corresponding to the X and Y samples, respectively, then the null is that F == G. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative F != G. In this case, the test statistic KS follows a two-sided Kolmogorov-Smirnov distribution. If ALT is ">", the one-sided alternative F > G is considered. Similarly for "<", the one-sided alternative F < G is considered. In this case, the test statistic KS has a one-sided Kolmogorov-Smirnov distribution. The default is the two-sided case. The p-value of the test is returned in PVAL. The third returned value, D, is the test statistic, the maximum vertical distance between the two cumulative distribution functions. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, K, DF] = kruskal_wallis_test (X1, ...) Perform a Kruskal-Wallis one-factor analysis of variance. Suppose a variable is observed for K > 1 different groups, and let X1, ..., XK be the corresponding data vectors. Under the null hypothesis that the ranks in the pooled sample are not affected by the group memberships, the test statistic K is approximately chi-square with DF = K - 1 degrees of freedom. If the data contains ties (some value appears more than once) K is divided by 1 - SUM_TIES / (N^3 - N) where SUM_TIES is the sum of T^2 - T over each group of ties where T is the number of ties in the group and N is the total number of values in the input data. For more info on this adjustment see William H. Kruskal and W. Allen Wallis, 'Use of Ranks in One-Criterion Variance Analysis', Journal of the American Statistical Association, Vol. 47, No. 260 (Dec 1952). The p-value (1 minus the CDF of this distribution at K) is returned in PVAL. If no output argument is given, the p-value is displayed. -- Function File: manova (X, G) Perform a one-way multivariate analysis of variance (MANOVA). The goal is to test whether the p-dimensional population means of data taken from K different groups are all equal. All data are assumed drawn independently from p-dimensional normal distributions with the same covariance matrix. The data matrix is given by X. As usual, rows are observations and columns are variables. The vector G specifies the corresponding group labels (e.g., numbers from 1 to K). The LR test statistic (Wilks' Lambda) and approximate p-values are computed and displayed. See also: *note anova: XREFanova. -- Function File: [PVAL, CHISQ, DF] = mcnemar_test (X) For a square contingency table X of data cross-classified on the row and column variables, McNemar's test can be used for testing the null hypothesis of symmetry of the classification probabilities. Under the null, CHISQ is approximately distributed as chisquare with DF degrees of freedom. The p-value (1 minus the CDF of this distribution at CHISQ) is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, Z] = prop_test_2 (X1, N1, X2, N2, ALT) If X1 and N1 are the counts of successes and trials in one sample, and X2 and N2 those in a second one, test the null hypothesis that the success probabilities P1 and P2 are the same. Under the null, the test statistic Z approximately follows a standard normal distribution. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative P1 != P2. If ALT is ">", the one-sided alternative P1 > P2 is used. Similarly for "<", the one-sided alternative P1 < P2 is used. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, CHISQ] = run_test (X) Perform a chi-square test with 6 degrees of freedom based on the upward runs in the columns of X. 'run_test' can be used to decide whether X contains independent data. The p-value of the test is returned in PVAL. If no output argument is given, the p-value is displayed. -- Function File: [PVAL, B, N] = sign_test (X, Y, ALT) For two matched-pair samples X and Y, perform a sign test of the null hypothesis PROB (X > Y) == PROB (X < Y) == 1/2. Under the null, the test statistic B roughly follows a binomial distribution with parameters 'N = sum (X != Y)' and P = 1/2. With the optional argument 'alt', the alternative of interest can be selected. If ALT is "!=" or "<>", the null hypothesis is tested against the two-sided alternative PROB (X < Y) != 1/2. If ALT is ">", the one-sided alternative PROB (X > Y) > 1/2 ("x is stochastically greater than y") is considered. Similarly for "<", the one-sided alternative PROB (X > Y) < 1/2 ("x is stochastically less than y") is considered. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, T, DF] = t_test (X, M, ALT) For a sample X from a normal distribution with unknown mean and variance, perform a t-test of the null hypothesis 'mean (X) == M'. Under the null, the test statistic T follows a Student distribution with 'DF = length (X) - 1' degrees of freedom. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'mean (X) != M'. If ALT is ">", the one-sided alternative 'mean (X) > M' is considered. Similarly for "<", the one-sided alternative 'mean (X) < M' is considered. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, T, DF] = t_test_2 (X, Y, ALT) For two samples x and y from normal distributions with unknown means and unknown equal variances, perform a two-sample t-test of the null hypothesis of equal means. Under the null, the test statistic T follows a Student distribution with DF degrees of freedom. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'mean (X) != mean (Y)'. If ALT is ">", the one-sided alternative 'mean (X) > mean (Y)' is used. Similarly for "<", the one-sided alternative 'mean (X) < mean (Y)' is used. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, T, DF] = t_test_regression (Y, X, RR, R, ALT) Perform a t test for the null hypothesis 'RR * B = R' in a classical normal regression model 'Y = X * B + E'. Under the null, the test statistic T follows a T distribution with DF degrees of freedom. If R is omitted, a value of 0 is assumed. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'RR * B != R'. If ALT is ">", the one-sided alternative 'RR * B > R' is used. Similarly for "<", the one-sided alternative 'RR * B < R' is used. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, Z] = u_test (X, Y, ALT) For two samples X and Y, perform a Mann-Whitney U-test of the null hypothesis PROB (X > Y) == 1/2 == PROB (X < Y). Under the null, the test statistic Z approximately follows a standard normal distribution. Note that this test is equivalent to the Wilcoxon rank-sum test. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative PROB (X > Y) != 1/2. If ALT is ">", the one-sided alternative PROB (X > Y) > 1/2 is considered. Similarly for "<", the one-sided alternative PROB (X > Y) < 1/2 is considered. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, F, DF_NUM, DF_DEN] = var_test (X, Y, ALT) For two samples X and Y from normal distributions with unknown means and unknown variances, perform an F-test of the null hypothesis of equal variances. Under the null, the test statistic F follows an F-distribution with DF_NUM and DF_DEN degrees of freedom. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'var (X) != var (Y)'. If ALT is ">", the one-sided alternative 'var (X) > var (Y)' is used. Similarly for "<", the one-sided alternative 'var (X) > var (Y)' is used. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, T, DF] = welch_test (X, Y, ALT) For two samples X and Y from normal distributions with unknown means and unknown and not necessarily equal variances, perform a Welch test of the null hypothesis of equal means. Under the null, the test statistic T approximately follows a Student distribution with DF degrees of freedom. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'mean (X) != M'. If ALT is ">", the one-sided alternative mean(x) > M is considered. Similarly for "<", the one-sided alternative mean(x) < M is considered. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, Z] = wilcoxon_test (X, Y, ALT) For two matched-pair sample vectors X and Y, perform a Wilcoxon signed-rank test of the null hypothesis PROB (X > Y) == 1/2. Under the null, the test statistic Z approximately follows a standard normal distribution when N > 25. *Caution:* This function assumes a normal distribution for Z and thus is invalid for N <= 25. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative PROB (X > Y) != 1/2. If alt is ">", the one-sided alternative PROB (X > Y) > 1/2 is considered. Similarly for "<", the one-sided alternative PROB (X > Y) < 1/2 is considered. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed. -- Function File: [PVAL, Z] = z_test (X, M, V, ALT) Perform a Z-test of the null hypothesis 'mean (X) == M' for a sample X from a normal distribution with unknown mean and known variance V. Under the null, the test statistic Z follows a standard normal distribution. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'mean (X) != M'. If ALT is ">", the one-sided alternative 'mean (X) > M' is considered. Similarly for "<", the one-sided alternative 'mean (X) < M' is considered. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed along with some information. -- Function File: [PVAL, Z] = z_test_2 (X, Y, V_X, V_Y, ALT) For two samples X and Y from normal distributions with unknown means and known variances V_X and V_Y, perform a Z-test of the hypothesis of equal means. Under the null, the test statistic Z follows a standard normal distribution. With the optional argument string ALT, the alternative of interest can be selected. If ALT is "!=" or "<>", the null is tested against the two-sided alternative 'mean (X) != mean (Y)'. If alt is ">", the one-sided alternative 'mean (X) > mean (Y)' is used. Similarly for "<", the one-sided alternative 'mean (X) < mean (Y)' is used. The default is the two-sided case. The p-value of the test is returned in PVAL. If no output argument is given, the p-value of the test is displayed along with some information.  File: octave.info, Node: Random Number Generation, Prev: Tests, Up: Statistics 26.7 Random Number Generation ============================= Octave can generate random numbers from a large number of distributions. The random number generators are based on the random number generators described in *note Special Utility Matrices::. The following table summarizes the available random number generators (in alphabetical order). Distribution Function ----------------------------------------------------- Beta Distribution 'betarnd' Binomial Distribution 'binornd' Cauchy Distribution 'cauchy_rnd' Chi-Square Distribution 'chi2rnd' Univariate Discrete 'discrete_rnd' Distribution Empirical Distribution 'empirical_rnd' Exponential Distribution 'exprnd' F Distribution 'frnd' Gamma Distribution 'gamrnd' Geometric Distribution 'geornd' Hypergeometric Distribution 'hygernd' Laplace Distribution 'laplace_rnd' Logistic Distribution 'logistic_rnd' Log-Normal Distribution 'lognrnd' Pascal Distribution 'nbinrnd' Univariate Normal 'normrnd' Distribution Poisson Distribution 'poissrnd' Standard Normal 'stdnormal_rnd' Distribution t (Student) Distribution 'trnd' Univariate Discrete 'unidrnd' Distribution Uniform Distribution 'unifrnd' Weibull Distribution 'wblrnd' Wiener Process 'wienrnd' -- Function File: betarnd (A, B) -- Function File: betarnd (A, B, R) -- Function File: betarnd (A, B, R, C, ...) -- Function File: betarnd (A, B, [SZ]) Return a matrix of random samples from the Beta distribution with parameters A and B. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of A and B. -- Function File: binornd (N, P) -- Function File: binornd (N, P, R) -- Function File: binornd (N, P, R, C, ...) -- Function File: binornd (N, P, [SZ]) Return a matrix of random samples from the binomial distribution with parameters N and P, where N is the number of trials and P is the probability of success. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of N and P. -- Function File: cauchy_rnd (LOCATION, SCALE) -- Function File: cauchy_rnd (LOCATION, SCALE, R) -- Function File: cauchy_rnd (LOCATION, SCALE, R, C, ...) -- Function File: cauchy_rnd (LOCATION, SCALE, [SZ]) Return a matrix of random samples from the Cauchy distribution with parameters LOCATION and SCALE. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of LOCATION and SCALE. -- Function File: chi2rnd (N) -- Function File: chi2rnd (N, R) -- Function File: chi2rnd (N, R, C, ...) -- Function File: chi2rnd (N, [SZ]) Return a matrix of random samples from the chi-square distribution with N degrees of freedom. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the size of N. -- Function File: discrete_rnd (V, P) -- Function File: discrete_rnd (V, P, R) -- Function File: discrete_rnd (V, P, R, C, ...) -- Function File: discrete_rnd (V, P, [SZ]) Return a matrix of random samples from the univariate distribution which assumes the values in V with probabilities P. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of V and P. -- Function File: empirical_rnd (DATA) -- Function File: empirical_rnd (DATA, R) -- Function File: empirical_rnd (DATA, R, C, ...) -- Function File: empirical_rnd (DATA, [SZ]) Return a matrix of random samples from the empirical distribution obtained from the univariate sample DATA. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is a random ordering of the sample DATA. -- Function File: exprnd (LAMBDA) -- Function File: exprnd (LAMBDA, R) -- Function File: exprnd (LAMBDA, R, C, ...) -- Function File: exprnd (LAMBDA, [SZ]) Return a matrix of random samples from the exponential distribution with mean LAMBDA. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the size of LAMBDA. -- Function File: frnd (M, N) -- Function File: frnd (M, N, R) -- Function File: frnd (M, N, R, C, ...) -- Function File: frnd (M, N, [SZ]) Return a matrix of random samples from the F distribution with M and N degrees of freedom. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of M and N. -- Function File: gamrnd (A, B) -- Function File: gamrnd (A, B, R) -- Function File: gamrnd (A, B, R, C, ...) -- Function File: gamrnd (A, B, [SZ]) Return a matrix of random samples from the Gamma distribution with shape parameter A and scale B. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of A and B. -- Function File: geornd (P) -- Function File: geornd (P, R) -- Function File: geornd (P, R, C, ...) -- Function File: geornd (P, [SZ]) Return a matrix of random samples from the geometric distribution with parameter P. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the size of P. The geometric distribution models the number of failures (X-1) of a Bernoulli trial with probability P before the first success (X). -- Function File: hygernd (T, M, N) -- Function File: hygernd (T, M, N, R) -- Function File: hygernd (T, M, N, R, C, ...) -- Function File: hygernd (T, M, N, [SZ]) Return a matrix of random samples from the hypergeometric distribution with parameters T, M, and N. The parameters T, M, and N must be positive integers with M and N not greater than T. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of T, M, and N. -- Function File: laplace_rnd (R) -- Function File: laplace_rnd (R, C, ...) -- Function File: laplace_rnd ([SZ]) Return a matrix of random samples from the Laplace distribution. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. -- Function File: logistic_rnd (R) -- Function File: logistic_rnd (R, C, ...) -- Function File: logistic_rnd ([SZ]) Return a matrix of random samples from the logistic distribution. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. -- Function File: lognrnd (MU, SIGMA) -- Function File: lognrnd (MU, SIGMA, R) -- Function File: lognrnd (MU, SIGMA, R, C, ...) -- Function File: lognrnd (MU, SIGMA, [SZ]) Return a matrix of random samples from the lognormal distribution with parameters MU and SIGMA. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of MU and SIGMA. -- Function File: nbinrnd (N, P) -- Function File: nbinrnd (N, P, R) -- Function File: nbinrnd (N, P, R, C, ...) -- Function File: nbinrnd (N, P, [SZ]) Return a matrix of random samples from the negative binomial distribution with parameters N and P. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of N and P. -- Function File: normrnd (MU, SIGMA) -- Function File: normrnd (MU, SIGMA, R) -- Function File: normrnd (MU, SIGMA, R, C, ...) -- Function File: normrnd (MU, SIGMA, [SZ]) Return a matrix of random samples from the normal distribution with parameters mean MU and standard deviation SIGMA. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of MU and SIGMA. -- Function File: poissrnd (LAMBDA) -- Function File: poissrnd (LAMBDA, R) -- Function File: poissrnd (LAMBDA, R, C, ...) -- Function File: poissrnd (LAMBDA, [SZ]) Return a matrix of random samples from the Poisson distribution with parameter LAMBDA. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the size of LAMBDA. -- Function File: stdnormal_rnd (R) -- Function File: stdnormal_rnd (R, C, ...) -- Function File: stdnormal_rnd ([SZ]) Return a matrix of random samples from the standard normal distribution (mean = 0, standard deviation = 1). When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. -- Function File: trnd (N) -- Function File: trnd (N, R) -- Function File: trnd (N, R, C, ...) -- Function File: trnd (N, [SZ]) Return a matrix of random samples from the t (Student) distribution with N degrees of freedom. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the size of N. -- Function File: unidrnd (N) -- Function File: unidrnd (N, R) -- Function File: unidrnd (N, R, C, ...) -- Function File: unidrnd (N, [SZ]) Return a matrix of random samples from the discrete uniform distribution which assumes the integer values 1-N with equal probability. N may be a scalar or a multi-dimensional array. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the size of N. -- Function File: unifrnd (A, B) -- Function File: unifrnd (A, B, R) -- Function File: unifrnd (A, B, R, C, ...) -- Function File: unifrnd (A, B, [SZ]) Return a matrix of random samples from the uniform distribution on [A, B]. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of A and B. -- Function File: wblrnd (SCALE, SHAPE) -- Function File: wblrnd (SCALE, SHAPE, R) -- Function File: wblrnd (SCALE, SHAPE, R, C, ...) -- Function File: wblrnd (SCALE, SHAPE, [SZ]) Return a matrix of random samples from the Weibull distribution with parameters SCALE and SHAPE. When called with a single size argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The size may also be specified with a vector of dimensions SZ. If no size arguments are given then the result matrix is the common size of SCALE and SHAPE. -- Function File: wienrnd (T, D, N) Return a simulated realization of the D-dimensional Wiener Process on the interval [0, T]. If D is omitted, D = 1 is used. The first column of the return matrix contains time, the remaining columns contain the Wiener process. The optional parameter N defines the number of summands used for simulating the process over an interval of length 1. If N is omitted, N = 1000 is used.  File: octave.info, Node: Sets, Next: Polynomial Manipulations, Prev: Statistics, Up: Top 27 Sets ******* Octave has a number of functions for managing sets of data. A set is defined as a collection of unique elements and is typically represented by a vector of numbers sorted in ascending order. Any vector or matrix can be converted to a set by removing duplicates through the use of the 'unique' function. However, it isn't necessary to explicitly create a set as all of the functions which operate on sets will convert their input to a set before proceeding. -- Function File: unique (X) -- Function File: unique (X, "rows") -- Function File: [Y, I, J] = unique (...) -- Function File: [Y, I, J] = unique (..., "first") -- Function File: [Y, I, J] = unique (..., "last") Return the unique elements of X sorted in ascending order. If the input X is a column vector then return a column vector; Otherwise, return a row vector. X may also be a cell array of strings. If the optional argument "rows" is given then return the unique rows of X sorted in ascending order. The input must be a 2-D matrix to use this option. If requested, return index vectors I and J such that 'Y = X(I)' and 'X = Y(J)'. Additionally, if I is a requested output then one of "first" or "last" may be given as an input. If "last" is specified, return the highest possible indices in I, otherwise, if "first" is specified, return the lowest. The default is "last". See also: *note union: XREFunion, *note intersect: XREFintersect, *note setdiff: XREFsetdiff, *note setxor: XREFsetxor, *note ismember: XREFismember. * Menu: * Set Operations::