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
18package org.apache.commons.math.linear;
19
20import org.apache.commons.math.util.FastMath;
21
22
23/**
24 * Class transforming any matrix to bi-diagonal shape.
25 * <p>Any m &times; n matrix A can be written as the product of three matrices:
26 * A = U &times; B &times; V<sup>T</sup> with U an m &times; m orthogonal matrix,
27 * B an m &times; n bi-diagonal matrix (lower diagonal if m &lt; n, upper diagonal
28 * otherwise), and V an n &times; n orthogonal matrix.</p>
29 * <p>Transformation to bi-diagonal shape is often not a goal by itself, but it is
30 * an intermediate step in more general decomposition algorithms like {@link
31 * SingularValueDecomposition Singular Value Decomposition}. This class is therefore
32 * intended for internal use by the library and is not public. As a consequence of
33 * this explicitly limited scope, many methods directly returns references to
34 * internal arrays, not copies.</p>
35 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
36 * @since 2.0
37 */
38class BiDiagonalTransformer {
39
40    /** Householder vectors. */
41    private final double householderVectors[][];
42
43    /** Main diagonal. */
44    private final double[] main;
45
46    /** Secondary diagonal. */
47    private final double[] secondary;
48
49    /** Cached value of U. */
50    private RealMatrix cachedU;
51
52    /** Cached value of B. */
53    private RealMatrix cachedB;
54
55    /** Cached value of V. */
56    private RealMatrix cachedV;
57
58    /**
59     * Build the transformation to bi-diagonal shape of a matrix.
60     * @param matrix the matrix to transform.
61     */
62    public BiDiagonalTransformer(RealMatrix matrix) {
63
64        final int m = matrix.getRowDimension();
65        final int n = matrix.getColumnDimension();
66        final int p = FastMath.min(m, n);
67        householderVectors = matrix.getData();
68        main      = new double[p];
69        secondary = new double[p - 1];
70        cachedU   = null;
71        cachedB   = null;
72        cachedV   = null;
73
74        // transform matrix
75        if (m >= n) {
76            transformToUpperBiDiagonal();
77        } else {
78            transformToLowerBiDiagonal();
79        }
80
81    }
82
83    /**
84     * Returns the matrix U of the transform.
85     * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
86     * @return the U matrix
87     */
88    public RealMatrix getU() {
89
90        if (cachedU == null) {
91
92            final int m = householderVectors.length;
93            final int n = householderVectors[0].length;
94            final int p = main.length;
95            final int diagOffset    = (m >= n) ? 0 : 1;
96            final double[] diagonal = (m >= n) ? main : secondary;
97            cachedU = MatrixUtils.createRealMatrix(m, m);
98
99            // fill up the part of the matrix not affected by Householder transforms
100            for (int k = m - 1; k >= p; --k) {
101                cachedU.setEntry(k, k, 1);
102            }
103
104            // build up first part of the matrix by applying Householder transforms
105            for (int k = p - 1; k >= diagOffset; --k) {
106                final double[] hK = householderVectors[k];
107                cachedU.setEntry(k, k, 1);
108                if (hK[k - diagOffset] != 0.0) {
109                    for (int j = k; j < m; ++j) {
110                        double alpha = 0;
111                        for (int i = k; i < m; ++i) {
112                            alpha -= cachedU.getEntry(i, j) * householderVectors[i][k - diagOffset];
113                        }
114                        alpha /= diagonal[k - diagOffset] * hK[k - diagOffset];
115
116                        for (int i = k; i < m; ++i) {
117                            cachedU.addToEntry(i, j, -alpha * householderVectors[i][k - diagOffset]);
118                        }
119                    }
120                }
121            }
122            if (diagOffset > 0) {
123                cachedU.setEntry(0, 0, 1);
124            }
125
126        }
127
128        // return the cached matrix
129        return cachedU;
130
131    }
132
133    /**
134     * Returns the bi-diagonal matrix B of the transform.
135     * @return the B matrix
136     */
137    public RealMatrix getB() {
138
139        if (cachedB == null) {
140
141            final int m = householderVectors.length;
142            final int n = householderVectors[0].length;
143            cachedB = MatrixUtils.createRealMatrix(m, n);
144            for (int i = 0; i < main.length; ++i) {
145                cachedB.setEntry(i, i, main[i]);
146                if (m < n) {
147                    if (i > 0) {
148                        cachedB.setEntry(i, i - 1, secondary[i - 1]);
149                    }
150                } else {
151                    if (i < main.length - 1) {
152                        cachedB.setEntry(i, i + 1, secondary[i]);
153                    }
154                }
155            }
156
157        }
158
159        // return the cached matrix
160        return cachedB;
161
162    }
163
164    /**
165     * Returns the matrix V of the transform.
166     * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
167     * @return the V matrix
168     */
169    public RealMatrix getV() {
170
171        if (cachedV == null) {
172
173            final int m = householderVectors.length;
174            final int n = householderVectors[0].length;
175            final int p = main.length;
176            final int diagOffset    = (m >= n) ? 1 : 0;
177            final double[] diagonal = (m >= n) ? secondary : main;
178            cachedV = MatrixUtils.createRealMatrix(n, n);
179
180            // fill up the part of the matrix not affected by Householder transforms
181            for (int k = n - 1; k >= p; --k) {
182                cachedV.setEntry(k, k, 1);
183            }
184
185            // build up first part of the matrix by applying Householder transforms
186            for (int k = p - 1; k >= diagOffset; --k) {
187                final double[] hK = householderVectors[k - diagOffset];
188                cachedV.setEntry(k, k, 1);
189                if (hK[k] != 0.0) {
190                    for (int j = k; j < n; ++j) {
191                        double beta = 0;
192                        for (int i = k; i < n; ++i) {
193                            beta -= cachedV.getEntry(i, j) * hK[i];
194                        }
195                        beta /= diagonal[k - diagOffset] * hK[k];
196
197                        for (int i = k; i < n; ++i) {
198                            cachedV.addToEntry(i, j, -beta * hK[i]);
199                        }
200                    }
201                }
202            }
203            if (diagOffset > 0) {
204                cachedV.setEntry(0, 0, 1);
205            }
206
207        }
208
209        // return the cached matrix
210        return cachedV;
211
212    }
213
214    /**
215     * Get the Householder vectors of the transform.
216     * <p>Note that since this class is only intended for internal use,
217     * it returns directly a reference to its internal arrays, not a copy.</p>
218     * @return the main diagonal elements of the B matrix
219     */
220    double[][] getHouseholderVectorsRef() {
221        return householderVectors;
222    }
223
224    /**
225     * Get the main diagonal elements of the matrix B of the transform.
226     * <p>Note that since this class is only intended for internal use,
227     * it returns directly a reference to its internal arrays, not a copy.</p>
228     * @return the main diagonal elements of the B matrix
229     */
230    double[] getMainDiagonalRef() {
231        return main;
232    }
233
234    /**
235     * Get the secondary diagonal elements of the matrix B of the transform.
236     * <p>Note that since this class is only intended for internal use,
237     * it returns directly a reference to its internal arrays, not a copy.</p>
238     * @return the secondary diagonal elements of the B matrix
239     */
240    double[] getSecondaryDiagonalRef() {
241        return secondary;
242    }
243
244    /**
245     * Check if the matrix is transformed to upper bi-diagonal.
246     * @return true if the matrix is transformed to upper bi-diagonal
247     */
248    boolean isUpperBiDiagonal() {
249        return householderVectors.length >=  householderVectors[0].length;
250    }
251
252    /**
253     * Transform original matrix to upper bi-diagonal form.
254     * <p>Transformation is done using alternate Householder transforms
255     * on columns and rows.</p>
256     */
257    private void transformToUpperBiDiagonal() {
258
259        final int m = householderVectors.length;
260        final int n = householderVectors[0].length;
261        for (int k = 0; k < n; k++) {
262
263            //zero-out a column
264            double xNormSqr = 0;
265            for (int i = k; i < m; ++i) {
266                final double c = householderVectors[i][k];
267                xNormSqr += c * c;
268            }
269            final double[] hK = householderVectors[k];
270            final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
271            main[k] = a;
272            if (a != 0.0) {
273                hK[k] -= a;
274                for (int j = k + 1; j < n; ++j) {
275                    double alpha = 0;
276                    for (int i = k; i < m; ++i) {
277                        final double[] hI = householderVectors[i];
278                        alpha -= hI[j] * hI[k];
279                    }
280                    alpha /= a * householderVectors[k][k];
281                    for (int i = k; i < m; ++i) {
282                        final double[] hI = householderVectors[i];
283                        hI[j] -= alpha * hI[k];
284                    }
285                }
286            }
287
288            if (k < n - 1) {
289                //zero-out a row
290                xNormSqr = 0;
291                for (int j = k + 1; j < n; ++j) {
292                    final double c = hK[j];
293                    xNormSqr += c * c;
294                }
295                final double b = (hK[k + 1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
296                secondary[k] = b;
297                if (b != 0.0) {
298                    hK[k + 1] -= b;
299                    for (int i = k + 1; i < m; ++i) {
300                        final double[] hI = householderVectors[i];
301                        double beta = 0;
302                        for (int j = k + 1; j < n; ++j) {
303                            beta -= hI[j] * hK[j];
304                        }
305                        beta /= b * hK[k + 1];
306                        for (int j = k + 1; j < n; ++j) {
307                            hI[j] -= beta * hK[j];
308                        }
309                    }
310                }
311            }
312
313        }
314    }
315
316    /**
317     * Transform original matrix to lower bi-diagonal form.
318     * <p>Transformation is done using alternate Householder transforms
319     * on rows and columns.</p>
320     */
321    private void transformToLowerBiDiagonal() {
322
323        final int m = householderVectors.length;
324        final int n = householderVectors[0].length;
325        for (int k = 0; k < m; k++) {
326
327            //zero-out a row
328            final double[] hK = householderVectors[k];
329            double xNormSqr = 0;
330            for (int j = k; j < n; ++j) {
331                final double c = hK[j];
332                xNormSqr += c * c;
333            }
334            final double a = (hK[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
335            main[k] = a;
336            if (a != 0.0) {
337                hK[k] -= a;
338                for (int i = k + 1; i < m; ++i) {
339                    final double[] hI = householderVectors[i];
340                    double alpha = 0;
341                    for (int j = k; j < n; ++j) {
342                        alpha -= hI[j] * hK[j];
343                    }
344                    alpha /= a * householderVectors[k][k];
345                    for (int j = k; j < n; ++j) {
346                        hI[j] -= alpha * hK[j];
347                    }
348                }
349            }
350
351            if (k < m - 1) {
352                //zero-out a column
353                final double[] hKp1 = householderVectors[k + 1];
354                xNormSqr = 0;
355                for (int i = k + 1; i < m; ++i) {
356                    final double c = householderVectors[i][k];
357                    xNormSqr += c * c;
358                }
359                final double b = (hKp1[k] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
360                secondary[k] = b;
361                if (b != 0.0) {
362                    hKp1[k] -= b;
363                    for (int j = k + 1; j < n; ++j) {
364                        double beta = 0;
365                        for (int i = k + 1; i < m; ++i) {
366                            final double[] hI = householderVectors[i];
367                            beta -= hI[j] * hI[k];
368                        }
369                        beta /= b * hKp1[k];
370                        for (int i = k + 1; i < m; ++i) {
371                            final double[] hI = householderVectors[i];
372                            hI[j] -= beta * hI[k];
373                        }
374                    }
375                }
376            }
377
378        }
379    }
380
381}
382