1 | /*
|
---|
2 | * Licensed to the Apache Software Foundation (ASF) under one or more
|
---|
3 | * contributor license agreements. See the NOTICE file distributed with
|
---|
4 | * this work for additional information regarding copyright ownership.
|
---|
5 | * The ASF licenses this file to You under the Apache License, Version 2.0
|
---|
6 | * (the "License"); you may not use this file except in compliance with
|
---|
7 | * the License. You may obtain a copy of the License at
|
---|
8 | *
|
---|
9 | * http://www.apache.org/licenses/LICENSE-2.0
|
---|
10 | *
|
---|
11 | * Unless required by applicable law or agreed to in writing, software
|
---|
12 | * distributed under the License is distributed on an "AS IS" BASIS,
|
---|
13 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
---|
14 | * See the License for the specific language governing permissions and
|
---|
15 | * limitations under the License.
|
---|
16 | */
|
---|
17 |
|
---|
18 | package agents.anac.y2019.harddealer.math3.linear;
|
---|
19 |
|
---|
20 | import agents.anac.y2019.harddealer.math3.util.FastMath;
|
---|
21 |
|
---|
22 | /**
|
---|
23 | * Calculates the rectangular Cholesky decomposition of a matrix.
|
---|
24 | * <p>The rectangular Cholesky decomposition of a real symmetric positive
|
---|
25 | * semidefinite matrix A consists of a rectangular matrix B with the same
|
---|
26 | * number of rows such that: A is almost equal to BB<sup>T</sup>, depending
|
---|
27 | * on a user-defined tolerance. In a sense, this is the square root of A.</p>
|
---|
28 | * <p>The difference with respect to the regular {@link CholeskyDecomposition}
|
---|
29 | * is that rows/columns may be permuted (hence the rectangular shape instead
|
---|
30 | * of the traditional triangular shape) and there is a threshold to ignore
|
---|
31 | * small diagonal elements. This is used for example to generate {@link
|
---|
32 | * agents.anac.y2019.harddealer.math3.random.CorrelatedRandomVectorGenerator correlated
|
---|
33 | * random n-dimensions vectors} in a p-dimension subspace (p < n).
|
---|
34 | * In other words, it allows generating random vectors from a covariance
|
---|
35 | * matrix that is only positive semidefinite, and not positive definite.</p>
|
---|
36 | * <p>Rectangular Cholesky decomposition is <em>not</em> suited for solving
|
---|
37 | * linear systems, so it does not provide any {@link DecompositionSolver
|
---|
38 | * decomposition solver}.</p>
|
---|
39 | *
|
---|
40 | * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
|
---|
41 | * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
|
---|
42 | * @since 2.0 (changed to concrete class in 3.0)
|
---|
43 | */
|
---|
44 | public class RectangularCholeskyDecomposition {
|
---|
45 |
|
---|
46 | /** Permutated Cholesky root of the symmetric positive semidefinite matrix. */
|
---|
47 | private final RealMatrix root;
|
---|
48 |
|
---|
49 | /** Rank of the symmetric positive semidefinite matrix. */
|
---|
50 | private int rank;
|
---|
51 |
|
---|
52 | /**
|
---|
53 | * Decompose a symmetric positive semidefinite matrix.
|
---|
54 | * <p>
|
---|
55 | * <b>Note:</b> this constructor follows the linpack method to detect dependent
|
---|
56 | * columns by proceeding with the Cholesky algorithm until a nonpositive diagonal
|
---|
57 | * element is encountered.
|
---|
58 | *
|
---|
59 | * @see <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf">
|
---|
60 | * Analysis of the Cholesky Decomposition of a Semi-definite Matrix</a>
|
---|
61 | *
|
---|
62 | * @param matrix Symmetric positive semidefinite matrix.
|
---|
63 | * @exception NonPositiveDefiniteMatrixException if the matrix is not
|
---|
64 | * positive semidefinite.
|
---|
65 | * @since 3.1
|
---|
66 | */
|
---|
67 | public RectangularCholeskyDecomposition(RealMatrix matrix)
|
---|
68 | throws NonPositiveDefiniteMatrixException {
|
---|
69 | this(matrix, 0);
|
---|
70 | }
|
---|
71 |
|
---|
72 | /**
|
---|
73 | * Decompose a symmetric positive semidefinite matrix.
|
---|
74 | *
|
---|
75 | * @param matrix Symmetric positive semidefinite matrix.
|
---|
76 | * @param small Diagonal elements threshold under which columns are
|
---|
77 | * considered to be dependent on previous ones and are discarded.
|
---|
78 | * @exception NonPositiveDefiniteMatrixException if the matrix is not
|
---|
79 | * positive semidefinite.
|
---|
80 | */
|
---|
81 | public RectangularCholeskyDecomposition(RealMatrix matrix, double small)
|
---|
82 | throws NonPositiveDefiniteMatrixException {
|
---|
83 |
|
---|
84 | final int order = matrix.getRowDimension();
|
---|
85 | final double[][] c = matrix.getData();
|
---|
86 | final double[][] b = new double[order][order];
|
---|
87 |
|
---|
88 | int[] index = new int[order];
|
---|
89 | for (int i = 0; i < order; ++i) {
|
---|
90 | index[i] = i;
|
---|
91 | }
|
---|
92 |
|
---|
93 | int r = 0;
|
---|
94 | for (boolean loop = true; loop;) {
|
---|
95 |
|
---|
96 | // find maximal diagonal element
|
---|
97 | int swapR = r;
|
---|
98 | for (int i = r + 1; i < order; ++i) {
|
---|
99 | int ii = index[i];
|
---|
100 | int isr = index[swapR];
|
---|
101 | if (c[ii][ii] > c[isr][isr]) {
|
---|
102 | swapR = i;
|
---|
103 | }
|
---|
104 | }
|
---|
105 |
|
---|
106 |
|
---|
107 | // swap elements
|
---|
108 | if (swapR != r) {
|
---|
109 | final int tmpIndex = index[r];
|
---|
110 | index[r] = index[swapR];
|
---|
111 | index[swapR] = tmpIndex;
|
---|
112 | final double[] tmpRow = b[r];
|
---|
113 | b[r] = b[swapR];
|
---|
114 | b[swapR] = tmpRow;
|
---|
115 | }
|
---|
116 |
|
---|
117 | // check diagonal element
|
---|
118 | int ir = index[r];
|
---|
119 | if (c[ir][ir] <= small) {
|
---|
120 |
|
---|
121 | if (r == 0) {
|
---|
122 | throw new NonPositiveDefiniteMatrixException(c[ir][ir], ir, small);
|
---|
123 | }
|
---|
124 |
|
---|
125 | // check remaining diagonal elements
|
---|
126 | for (int i = r; i < order; ++i) {
|
---|
127 | if (c[index[i]][index[i]] < -small) {
|
---|
128 | // there is at least one sufficiently negative diagonal element,
|
---|
129 | // the symmetric positive semidefinite matrix is wrong
|
---|
130 | throw new NonPositiveDefiniteMatrixException(c[index[i]][index[i]], i, small);
|
---|
131 | }
|
---|
132 | }
|
---|
133 |
|
---|
134 | // all remaining diagonal elements are close to zero, we consider we have
|
---|
135 | // found the rank of the symmetric positive semidefinite matrix
|
---|
136 | loop = false;
|
---|
137 |
|
---|
138 | } else {
|
---|
139 |
|
---|
140 | // transform the matrix
|
---|
141 | final double sqrt = FastMath.sqrt(c[ir][ir]);
|
---|
142 | b[r][r] = sqrt;
|
---|
143 | final double inverse = 1 / sqrt;
|
---|
144 | final double inverse2 = 1 / c[ir][ir];
|
---|
145 | for (int i = r + 1; i < order; ++i) {
|
---|
146 | final int ii = index[i];
|
---|
147 | final double e = inverse * c[ii][ir];
|
---|
148 | b[i][r] = e;
|
---|
149 | c[ii][ii] -= c[ii][ir] * c[ii][ir] * inverse2;
|
---|
150 | for (int j = r + 1; j < i; ++j) {
|
---|
151 | final int ij = index[j];
|
---|
152 | final double f = c[ii][ij] - e * b[j][r];
|
---|
153 | c[ii][ij] = f;
|
---|
154 | c[ij][ii] = f;
|
---|
155 | }
|
---|
156 | }
|
---|
157 |
|
---|
158 | // prepare next iteration
|
---|
159 | loop = ++r < order;
|
---|
160 | }
|
---|
161 | }
|
---|
162 |
|
---|
163 | // build the root matrix
|
---|
164 | rank = r;
|
---|
165 | root = MatrixUtils.createRealMatrix(order, r);
|
---|
166 | for (int i = 0; i < order; ++i) {
|
---|
167 | for (int j = 0; j < r; ++j) {
|
---|
168 | root.setEntry(index[i], j, b[i][j]);
|
---|
169 | }
|
---|
170 | }
|
---|
171 |
|
---|
172 | }
|
---|
173 |
|
---|
174 | /** Get the root of the covariance matrix.
|
---|
175 | * The root is the rectangular matrix <code>B</code> such that
|
---|
176 | * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
|
---|
177 | * @return root of the square matrix
|
---|
178 | * @see #getRank()
|
---|
179 | */
|
---|
180 | public RealMatrix getRootMatrix() {
|
---|
181 | return root;
|
---|
182 | }
|
---|
183 |
|
---|
184 | /** Get the rank of the symmetric positive semidefinite matrix.
|
---|
185 | * The r is the number of independent rows in the symmetric positive semidefinite
|
---|
186 | * matrix, it is also the number of columns of the rectangular
|
---|
187 | * matrix of the decomposition.
|
---|
188 | * @return r of the square matrix.
|
---|
189 | * @see #getRootMatrix()
|
---|
190 | */
|
---|
191 | public int getRank() {
|
---|
192 | return rank;
|
---|
193 | }
|
---|
194 |
|
---|
195 | }
|
---|