001/** 002 * Copyright (c) 2011, The University of Southampton and the individual contributors. 003 * All rights reserved. 004 * 005 * Redistribution and use in source and binary forms, with or without modification, 006 * are permitted provided that the following conditions are met: 007 * 008 * * Redistributions of source code must retain the above copyright notice, 009 * this list of conditions and the following disclaimer. 010 * 011 * * Redistributions in binary form must reproduce the above copyright notice, 012 * this list of conditions and the following disclaimer in the documentation 013 * and/or other materials provided with the distribution. 014 * 015 * * Neither the name of the University of Southampton nor the names of its 016 * contributors may be used to endorse or promote products derived from this 017 * software without specific prior written permission. 018 * 019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 029 */ 030package org.openimaj.math.matrix; 031 032import Jama.CholeskyDecomposition; 033import Jama.Matrix; 034 035/** 036 * Cholesky Update and Downdate 037 * 038 * @author Sina Samangooei (ss@ecs.soton.ac.uk) 039 */ 040public class UpdateableCholeskyDecomposition extends CholeskyDecomposition { 041 042 /** 043 * 044 */ 045 private static final long serialVersionUID = -7456377954967521480L; 046 047 /** 048 * @param m 049 */ 050 public UpdateableCholeskyDecomposition(Matrix m) { 051 super(m); 052 } 053 054 055 /** 056 * @param x 057 */ 058 public void choldowndate(double[] x) { 059 choldowndate(x,true); 060 } 061 062 /** 063 * see {@link UpdateableCholeskyDecomposition#choldowndate(double[][], double[])} 064 * @param x 065 * @param b 066 */ 067 public void choldowndate(double[] x, boolean b) { 068 if(b) x = x.clone(); 069 Matrix L = this.getL(); 070 // work is done on an upper triangular matrix 071 double[][] data = L.transpose().getArray(); 072 choldowndate(data, x); 073 // Make the output lower triangular again 074 int Ll = L.getRowDimension(); 075 L.setMatrix(0, Ll-1, 0, Ll-1, new Matrix(data, Ll, Ll).transpose()); 076 } 077 078 079 /** 080 * See {@link UpdateableCholeskyDecomposition#cholupdate(double[][], double[])} 081 * @param x 082 */ 083 public void cholupdate(double[] x){ 084 cholupdate(x,true); 085 } 086 087 /** 088 * See {@link #cholupdate(double[][], double[])} 089 * @param x 090 * @param copyX copy x, x is used as a workspace 091 */ 092 public void cholupdate(double[] x, boolean copyX){ 093 if(copyX){ 094 x = x.clone(); 095 } 096 Matrix L = this.getL(); 097 // work is done on an upper triangular matrix 098 double[][] data = L.transpose().getArray(); 099 cholupdate(data, x); 100 // Make the output lower triangular again 101 int Ll = L.getRowDimension(); 102 L.setMatrix(0, Ll-1, 0, Ll-1, new Matrix(data, Ll, Ll).transpose()); 103 } 104 /** 105 * Updates L such that the matrix A where: 106 * A = LL* 107 * becomes 108 * A = A + xx* 109 * @param Larr the upper triangular matrix to update 110 * @param x the vector to add 111 */ 112 public static void cholupdate(double[][] Larr, double[] x){ 113// x = x'; 114 for (int k = 0; k < x.length; k++) { 115 double Lkk = Larr[k][k]; 116 double xk = x[k]; 117 double r = Math.sqrt(Lkk*Lkk + xk*xk); 118 double c = r / Lkk; 119 double s = xk / Lkk; 120 Larr[k][k] = r; 121 updateL(k,Larr,x,s,c); 122 updateX(k,Larr,x,s,c); 123 } 124 } 125 126 /** 127 * Updates L such that the matrix A where: 128 * A = LL* 129 * becomes 130 * A = A - xx* 131 * @param Larr the upper triangular matrix to update 132 * @param x the vector to add 133 */ 134 public static void choldowndate(double[][] Larr, double[] x){ 135// x = x'; 136 for (int k = 0; k < x.length; k++) { 137 double Lkk = Larr[k][k]; 138 double xk = x[k]; 139 double r = Math.sqrt(Lkk*Lkk - xk*xk); 140 double c = r / Lkk; 141 double s = xk / Lkk; 142 Larr[k][k] = r; 143 downdateL(k,Larr,x,s,c); 144 updateX(k,Larr,x,s,c); 145 } 146 } 147 148 /** 149 * x(k+1:p) = c*x(k+1:p) - s*L(k, k+1:p); 150 * 151 * p = x.length 152 * 153 * @param k 154 * @param Larr 155 * @param x 156 * @param s 157 * @param c 158 */ 159 private static void updateX(int k, double[][] Larr, double[] x, double s, double c) { 160 for (int i = k+1; i < x.length; i++) { 161 x[i] = c*x[i] - s * Larr[k][i]; 162 } 163 } 164 165 /** 166 * L(k,k+1:p) = (L(k,k+1:p) + s*x(k+1:p)) / c; 167 * 168 * p = x.length 169 * @param k 170 * @param Larr 171 * @param s 172 * @param x 173 * @param c 174 */ 175 private static void updateL(int k, double[][] Larr, double[] x, double s, double c) { 176 for (int i = k+1; i < x.length; i++) { 177 Larr[k][i] = (Larr[k][i] + s * x[i])/c; 178 } 179 } 180 181 /** 182 * L(k,k+1:p) = (L(k,k+1:p) - s*x(k+1:p)) / c; 183 * 184 * p = x.length 185 * @param k 186 * @param Larr 187 * @param s 188 * @param x 189 * @param c 190 */ 191 private static void downdateL(int k, double[][] Larr, double[] x, double s, double c) { 192 for (int i = k+1; i < x.length; i++) { 193 Larr[k][i] = (Larr[k][i] - s * x[i])/c; 194 } 195 } 196}