Mercurial > hg > de.mpg.mpiwg.itgroup.digilib.plugin

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/linear/SingularValueDecompositionImpl.html @ 13:cbf34dd4d7e6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
commons-math-2.1 added
author dwinter
date Tue, 04 Jan 2011 10:02:07 +0100
parents
children
comparison
equal deleted inserted replaced
12:970d26a94fb7 13:cbf34dd4d7e6
1 <HTML>
2 <BODY BGCOLOR="white">
3 <PRE>
4 <FONT color="green">001</FONT> /*<a name="line.1"></a>
5 <FONT color="green">002</FONT> * Licensed to the Apache Software Foundation (ASF) under one or more<a name="line.2"></a>
6 <FONT color="green">003</FONT> * contributor license agreements. See the NOTICE file distributed with<a name="line.3"></a>
7 <FONT color="green">004</FONT> * this work for additional information regarding copyright ownership.<a name="line.4"></a>
8 <FONT color="green">005</FONT> * The ASF licenses this file to You under the Apache License, Version 2.0<a name="line.5"></a>
9 <FONT color="green">006</FONT> * (the "License"); you may not use this file except in compliance with<a name="line.6"></a>
10 <FONT color="green">007</FONT> * the License. You may obtain a copy of the License at<a name="line.7"></a>
11 <FONT color="green">008</FONT> *<a name="line.8"></a>
12 <FONT color="green">009</FONT> * http://www.apache.org/licenses/LICENSE-2.0<a name="line.9"></a>
13 <FONT color="green">010</FONT> *<a name="line.10"></a>
14 <FONT color="green">011</FONT> * Unless required by applicable law or agreed to in writing, software<a name="line.11"></a>
15 <FONT color="green">012</FONT> * distributed under the License is distributed on an "AS IS" BASIS,<a name="line.12"></a>
16 <FONT color="green">013</FONT> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.<a name="line.13"></a>
17 <FONT color="green">014</FONT> * See the License for the specific language governing permissions and<a name="line.14"></a>
18 <FONT color="green">015</FONT> * limitations under the License.<a name="line.15"></a>
19 <FONT color="green">016</FONT> */<a name="line.16"></a>
20 <FONT color="green">017</FONT> <a name="line.17"></a>
21 <FONT color="green">018</FONT> package org.apache.commons.math.linear;<a name="line.18"></a>
22 <FONT color="green">019</FONT> <a name="line.19"></a>
23 <FONT color="green">020</FONT> import org.apache.commons.math.MathRuntimeException;<a name="line.20"></a>
24 <FONT color="green">021</FONT> <a name="line.21"></a>
25 <FONT color="green">022</FONT> /**<a name="line.22"></a>
26 <FONT color="green">023</FONT> * Calculates the compact Singular Value Decomposition of a matrix.<a name="line.23"></a>
27 <FONT color="green">024</FONT> * &lt;p&gt;<a name="line.24"></a>
28 <FONT color="green">025</FONT> * The Singular Value Decomposition of matrix A is a set of three matrices: U,<a name="line.25"></a>
29 <FONT color="green">026</FONT> * &amp;Sigma; and V such that A = U &amp;times; &amp;Sigma; &amp;times; V&lt;sup&gt;T&lt;/sup&gt;. Let A be<a name="line.26"></a>
30 <FONT color="green">027</FONT> * a m &amp;times; n matrix, then U is a m &amp;times; p orthogonal matrix, &amp;Sigma; is a<a name="line.27"></a>
31 <FONT color="green">028</FONT> * p &amp;times; p diagonal matrix with positive or null elements, V is a p &amp;times;<a name="line.28"></a>
32 <FONT color="green">029</FONT> * n orthogonal matrix (hence V&lt;sup&gt;T&lt;/sup&gt; is also orthogonal) where<a name="line.29"></a>
33 <FONT color="green">030</FONT> * p=min(m,n).<a name="line.30"></a>
34 <FONT color="green">031</FONT> * &lt;/p&gt;<a name="line.31"></a>
35 <FONT color="green">032</FONT> * @version $Revision: 912413 $ $Date: 2010-02-21 16:46:12 -0500 (Sun, 21 Feb 2010) $<a name="line.32"></a>
36 <FONT color="green">033</FONT> * @since 2.0<a name="line.33"></a>
37 <FONT color="green">034</FONT> */<a name="line.34"></a>
38 <FONT color="green">035</FONT> public class SingularValueDecompositionImpl implements<a name="line.35"></a>
39 <FONT color="green">036</FONT> SingularValueDecomposition {<a name="line.36"></a>
40 <FONT color="green">037</FONT> <a name="line.37"></a>
41 <FONT color="green">038</FONT> /** Number of rows of the initial matrix. */<a name="line.38"></a>
42 <FONT color="green">039</FONT> private int m;<a name="line.39"></a>
43 <FONT color="green">040</FONT> <a name="line.40"></a>
44 <FONT color="green">041</FONT> /** Number of columns of the initial matrix. */<a name="line.41"></a>
45 <FONT color="green">042</FONT> private int n;<a name="line.42"></a>
46 <FONT color="green">043</FONT> <a name="line.43"></a>
47 <FONT color="green">044</FONT> /** Eigen decomposition of the tridiagonal matrix. */<a name="line.44"></a>
48 <FONT color="green">045</FONT> private EigenDecomposition eigenDecomposition;<a name="line.45"></a>
49 <FONT color="green">046</FONT> <a name="line.46"></a>
50 <FONT color="green">047</FONT> /** Singular values. */<a name="line.47"></a>
51 <FONT color="green">048</FONT> private double[] singularValues;<a name="line.48"></a>
52 <FONT color="green">049</FONT> <a name="line.49"></a>
53 <FONT color="green">050</FONT> /** Cached value of U. */<a name="line.50"></a>
54 <FONT color="green">051</FONT> private RealMatrix cachedU;<a name="line.51"></a>
55 <FONT color="green">052</FONT> <a name="line.52"></a>
56 <FONT color="green">053</FONT> /** Cached value of U&lt;sup&gt;T&lt;/sup&gt;. */<a name="line.53"></a>
57 <FONT color="green">054</FONT> private RealMatrix cachedUt;<a name="line.54"></a>
58 <FONT color="green">055</FONT> <a name="line.55"></a>
59 <FONT color="green">056</FONT> /** Cached value of S. */<a name="line.56"></a>
60 <FONT color="green">057</FONT> private RealMatrix cachedS;<a name="line.57"></a>
61 <FONT color="green">058</FONT> <a name="line.58"></a>
62 <FONT color="green">059</FONT> /** Cached value of V. */<a name="line.59"></a>
63 <FONT color="green">060</FONT> private RealMatrix cachedV;<a name="line.60"></a>
64 <FONT color="green">061</FONT> <a name="line.61"></a>
65 <FONT color="green">062</FONT> /** Cached value of V&lt;sup&gt;T&lt;/sup&gt;. */<a name="line.62"></a>
66 <FONT color="green">063</FONT> private RealMatrix cachedVt;<a name="line.63"></a>
67 <FONT color="green">064</FONT> <a name="line.64"></a>
68 <FONT color="green">065</FONT> /**<a name="line.65"></a>
69 <FONT color="green">066</FONT> * Calculates the compact Singular Value Decomposition of the given matrix.<a name="line.66"></a>
70 <FONT color="green">067</FONT> * @param matrix<a name="line.67"></a>
71 <FONT color="green">068</FONT> * The matrix to decompose.<a name="line.68"></a>
72 <FONT color="green">069</FONT> * @exception InvalidMatrixException<a name="line.69"></a>
73 <FONT color="green">070</FONT> * (wrapping a<a name="line.70"></a>
74 <FONT color="green">071</FONT> * {@link org.apache.commons.math.ConvergenceException} if<a name="line.71"></a>
75 <FONT color="green">072</FONT> * algorithm fails to converge<a name="line.72"></a>
76 <FONT color="green">073</FONT> */<a name="line.73"></a>
77 <FONT color="green">074</FONT> public SingularValueDecompositionImpl(final RealMatrix matrix)<a name="line.74"></a>
78 <FONT color="green">075</FONT> throws InvalidMatrixException {<a name="line.75"></a>
79 <FONT color="green">076</FONT> <a name="line.76"></a>
80 <FONT color="green">077</FONT> m = matrix.getRowDimension();<a name="line.77"></a>
81 <FONT color="green">078</FONT> n = matrix.getColumnDimension();<a name="line.78"></a>
82 <FONT color="green">079</FONT> <a name="line.79"></a>
83 <FONT color="green">080</FONT> cachedU = null;<a name="line.80"></a>
84 <FONT color="green">081</FONT> cachedS = null;<a name="line.81"></a>
85 <FONT color="green">082</FONT> cachedV = null;<a name="line.82"></a>
86 <FONT color="green">083</FONT> cachedVt = null;<a name="line.83"></a>
87 <FONT color="green">084</FONT> <a name="line.84"></a>
88 <FONT color="green">085</FONT> double[][] localcopy = matrix.getData();<a name="line.85"></a>
89 <FONT color="green">086</FONT> double[][] matATA = new double[n][n];<a name="line.86"></a>
90 <FONT color="green">087</FONT> //<a name="line.87"></a>
91 <FONT color="green">088</FONT> // create A^T*A<a name="line.88"></a>
92 <FONT color="green">089</FONT> //<a name="line.89"></a>
93 <FONT color="green">090</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.90"></a>
94 <FONT color="green">091</FONT> for (int j = i; j &lt; n; j++) {<a name="line.91"></a>
95 <FONT color="green">092</FONT> matATA[i][j] = 0.0;<a name="line.92"></a>
96 <FONT color="green">093</FONT> for (int k = 0; k &lt; m; k++) {<a name="line.93"></a>
97 <FONT color="green">094</FONT> matATA[i][j] += localcopy[k][i] * localcopy[k][j];<a name="line.94"></a>
98 <FONT color="green">095</FONT> }<a name="line.95"></a>
99 <FONT color="green">096</FONT> matATA[j][i]=matATA[i][j];<a name="line.96"></a>
100 <FONT color="green">097</FONT> }<a name="line.97"></a>
101 <FONT color="green">098</FONT> }<a name="line.98"></a>
102 <FONT color="green">099</FONT> <a name="line.99"></a>
103 <FONT color="green">100</FONT> double[][] matAAT = new double[m][m];<a name="line.100"></a>
104 <FONT color="green">101</FONT> //<a name="line.101"></a>
105 <FONT color="green">102</FONT> // create A*A^T<a name="line.102"></a>
106 <FONT color="green">103</FONT> //<a name="line.103"></a>
107 <FONT color="green">104</FONT> for (int i = 0; i &lt; m; i++) {<a name="line.104"></a>
108 <FONT color="green">105</FONT> for (int j = i; j &lt; m; j++) {<a name="line.105"></a>
109 <FONT color="green">106</FONT> matAAT[i][j] = 0.0;<a name="line.106"></a>
110 <FONT color="green">107</FONT> for (int k = 0; k &lt; n; k++) {<a name="line.107"></a>
111 <FONT color="green">108</FONT> matAAT[i][j] += localcopy[i][k] * localcopy[j][k];<a name="line.108"></a>
112 <FONT color="green">109</FONT> }<a name="line.109"></a>
113 <FONT color="green">110</FONT> matAAT[j][i]=matAAT[i][j];<a name="line.110"></a>
114 <FONT color="green">111</FONT> }<a name="line.111"></a>
115 <FONT color="green">112</FONT> }<a name="line.112"></a>
116 <FONT color="green">113</FONT> int p;<a name="line.113"></a>
117 <FONT color="green">114</FONT> if (m&gt;=n) {<a name="line.114"></a>
118 <FONT color="green">115</FONT> p=n;<a name="line.115"></a>
119 <FONT color="green">116</FONT> // compute eigen decomposition of A^T*A<a name="line.116"></a>
120 <FONT color="green">117</FONT> eigenDecomposition = new EigenDecompositionImpl(<a name="line.117"></a>
121 <FONT color="green">118</FONT> new Array2DRowRealMatrix(matATA),1.0);<a name="line.118"></a>
122 <FONT color="green">119</FONT> singularValues = eigenDecomposition.getRealEigenvalues();<a name="line.119"></a>
123 <FONT color="green">120</FONT> cachedV = eigenDecomposition.getV();<a name="line.120"></a>
124 <FONT color="green">121</FONT> <a name="line.121"></a>
125 <FONT color="green">122</FONT> // compute eigen decomposition of A*A^T<a name="line.122"></a>
126 <FONT color="green">123</FONT> eigenDecomposition = new EigenDecompositionImpl(<a name="line.123"></a>
127 <FONT color="green">124</FONT> new Array2DRowRealMatrix(matAAT),1.0);<a name="line.124"></a>
128 <FONT color="green">125</FONT> cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);<a name="line.125"></a>
129 <FONT color="green">126</FONT> } else {<a name="line.126"></a>
130 <FONT color="green">127</FONT> p=m;<a name="line.127"></a>
131 <FONT color="green">128</FONT> // compute eigen decomposition of A*A^T<a name="line.128"></a>
132 <FONT color="green">129</FONT> eigenDecomposition = new EigenDecompositionImpl(<a name="line.129"></a>
133 <FONT color="green">130</FONT> new Array2DRowRealMatrix(matAAT),1.0);<a name="line.130"></a>
134 <FONT color="green">131</FONT> singularValues = eigenDecomposition.getRealEigenvalues();<a name="line.131"></a>
135 <FONT color="green">132</FONT> cachedU = eigenDecomposition.getV();<a name="line.132"></a>
136 <FONT color="green">133</FONT> <a name="line.133"></a>
137 <FONT color="green">134</FONT> // compute eigen decomposition of A^T*A<a name="line.134"></a>
138 <FONT color="green">135</FONT> eigenDecomposition = new EigenDecompositionImpl(<a name="line.135"></a>
139 <FONT color="green">136</FONT> new Array2DRowRealMatrix(matATA),1.0);<a name="line.136"></a>
140 <FONT color="green">137</FONT> cachedV = eigenDecomposition.getV().getSubMatrix(0,n-1,0,p-1);<a name="line.137"></a>
141 <FONT color="green">138</FONT> }<a name="line.138"></a>
142 <FONT color="green">139</FONT> for (int i = 0; i &lt; p; i++) {<a name="line.139"></a>
143 <FONT color="green">140</FONT> singularValues[i] = Math.sqrt(Math.abs(singularValues[i]));<a name="line.140"></a>
144 <FONT color="green">141</FONT> }<a name="line.141"></a>
145 <FONT color="green">142</FONT> // Up to this point, U and V are computed independently of each other.<a name="line.142"></a>
146 <FONT color="green">143</FONT> // There still an sign indetermination of each column of, say, U.<a name="line.143"></a>
147 <FONT color="green">144</FONT> // The sign is set such that A.V_i=sigma_i.U_i (i&lt;=p)<a name="line.144"></a>
148 <FONT color="green">145</FONT> // The right sign corresponds to a positive dot product of A.V_i and U_i<a name="line.145"></a>
149 <FONT color="green">146</FONT> for (int i = 0; i &lt; p; i++) {<a name="line.146"></a>
150 <FONT color="green">147</FONT> RealVector tmp = cachedU.getColumnVector(i);<a name="line.147"></a>
151 <FONT color="green">148</FONT> double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp);<a name="line.148"></a>
152 <FONT color="green">149</FONT> if (product&lt;0) {<a name="line.149"></a>
153 <FONT color="green">150</FONT> cachedU.setColumnVector(i, tmp.mapMultiply(-1.0));<a name="line.150"></a>
154 <FONT color="green">151</FONT> }<a name="line.151"></a>
155 <FONT color="green">152</FONT> }<a name="line.152"></a>
156 <FONT color="green">153</FONT> }<a name="line.153"></a>
157 <FONT color="green">154</FONT> <a name="line.154"></a>
158 <FONT color="green">155</FONT> /** {@inheritDoc} */<a name="line.155"></a>
159 <FONT color="green">156</FONT> public RealMatrix getU() throws InvalidMatrixException {<a name="line.156"></a>
160 <FONT color="green">157</FONT> // return the cached matrix<a name="line.157"></a>
161 <FONT color="green">158</FONT> return cachedU;<a name="line.158"></a>
162 <FONT color="green">159</FONT> <a name="line.159"></a>
163 <FONT color="green">160</FONT> }<a name="line.160"></a>
164 <FONT color="green">161</FONT> <a name="line.161"></a>
165 <FONT color="green">162</FONT> /** {@inheritDoc} */<a name="line.162"></a>
166 <FONT color="green">163</FONT> public RealMatrix getUT() throws InvalidMatrixException {<a name="line.163"></a>
167 <FONT color="green">164</FONT> <a name="line.164"></a>
168 <FONT color="green">165</FONT> if (cachedUt == null) {<a name="line.165"></a>
169 <FONT color="green">166</FONT> cachedUt = getU().transpose();<a name="line.166"></a>
170 <FONT color="green">167</FONT> }<a name="line.167"></a>
171 <FONT color="green">168</FONT> <a name="line.168"></a>
172 <FONT color="green">169</FONT> // return the cached matrix<a name="line.169"></a>
173 <FONT color="green">170</FONT> return cachedUt;<a name="line.170"></a>
174 <FONT color="green">171</FONT> <a name="line.171"></a>
175 <FONT color="green">172</FONT> }<a name="line.172"></a>
176 <FONT color="green">173</FONT> <a name="line.173"></a>
177 <FONT color="green">174</FONT> /** {@inheritDoc} */<a name="line.174"></a>
178 <FONT color="green">175</FONT> public RealMatrix getS() throws InvalidMatrixException {<a name="line.175"></a>
179 <FONT color="green">176</FONT> <a name="line.176"></a>
180 <FONT color="green">177</FONT> if (cachedS == null) {<a name="line.177"></a>
181 <FONT color="green">178</FONT> <a name="line.178"></a>
182 <FONT color="green">179</FONT> // cache the matrix for subsequent calls<a name="line.179"></a>
183 <FONT color="green">180</FONT> cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);<a name="line.180"></a>
184 <FONT color="green">181</FONT> <a name="line.181"></a>
185 <FONT color="green">182</FONT> }<a name="line.182"></a>
186 <FONT color="green">183</FONT> return cachedS;<a name="line.183"></a>
187 <FONT color="green">184</FONT> }<a name="line.184"></a>
188 <FONT color="green">185</FONT> <a name="line.185"></a>
189 <FONT color="green">186</FONT> /** {@inheritDoc} */<a name="line.186"></a>
190 <FONT color="green">187</FONT> public double[] getSingularValues() throws InvalidMatrixException {<a name="line.187"></a>
191 <FONT color="green">188</FONT> return singularValues.clone();<a name="line.188"></a>
192 <FONT color="green">189</FONT> }<a name="line.189"></a>
193 <FONT color="green">190</FONT> <a name="line.190"></a>
194 <FONT color="green">191</FONT> /** {@inheritDoc} */<a name="line.191"></a>
195 <FONT color="green">192</FONT> public RealMatrix getV() throws InvalidMatrixException {<a name="line.192"></a>
196 <FONT color="green">193</FONT> // return the cached matrix<a name="line.193"></a>
197 <FONT color="green">194</FONT> return cachedV;<a name="line.194"></a>
198 <FONT color="green">195</FONT> <a name="line.195"></a>
199 <FONT color="green">196</FONT> }<a name="line.196"></a>
200 <FONT color="green">197</FONT> <a name="line.197"></a>
201 <FONT color="green">198</FONT> /** {@inheritDoc} */<a name="line.198"></a>
202 <FONT color="green">199</FONT> public RealMatrix getVT() throws InvalidMatrixException {<a name="line.199"></a>
203 <FONT color="green">200</FONT> <a name="line.200"></a>
204 <FONT color="green">201</FONT> if (cachedVt == null) {<a name="line.201"></a>
205 <FONT color="green">202</FONT> cachedVt = getV().transpose();<a name="line.202"></a>
206 <FONT color="green">203</FONT> }<a name="line.203"></a>
207 <FONT color="green">204</FONT> <a name="line.204"></a>
208 <FONT color="green">205</FONT> // return the cached matrix<a name="line.205"></a>
209 <FONT color="green">206</FONT> return cachedVt;<a name="line.206"></a>
210 <FONT color="green">207</FONT> <a name="line.207"></a>
211 <FONT color="green">208</FONT> }<a name="line.208"></a>
212 <FONT color="green">209</FONT> <a name="line.209"></a>
213 <FONT color="green">210</FONT> /** {@inheritDoc} */<a name="line.210"></a>
214 <FONT color="green">211</FONT> public RealMatrix getCovariance(final double minSingularValue) {<a name="line.211"></a>
215 <FONT color="green">212</FONT> <a name="line.212"></a>
216 <FONT color="green">213</FONT> // get the number of singular values to consider<a name="line.213"></a>
217 <FONT color="green">214</FONT> final int p = singularValues.length;<a name="line.214"></a>
218 <FONT color="green">215</FONT> int dimension = 0;<a name="line.215"></a>
219 <FONT color="green">216</FONT> while ((dimension &lt; p) &amp;&amp; (singularValues[dimension] &gt;= minSingularValue)) {<a name="line.216"></a>
220 <FONT color="green">217</FONT> ++dimension;<a name="line.217"></a>
221 <FONT color="green">218</FONT> }<a name="line.218"></a>
222 <FONT color="green">219</FONT> <a name="line.219"></a>
223 <FONT color="green">220</FONT> if (dimension == 0) {<a name="line.220"></a>
224 <FONT color="green">221</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.221"></a>
225 <FONT color="green">222</FONT> "cutoff singular value is {0}, should be at most {1}",<a name="line.222"></a>
226 <FONT color="green">223</FONT> minSingularValue, singularValues[0]);<a name="line.223"></a>
227 <FONT color="green">224</FONT> }<a name="line.224"></a>
228 <FONT color="green">225</FONT> <a name="line.225"></a>
229 <FONT color="green">226</FONT> final double[][] data = new double[dimension][p];<a name="line.226"></a>
230 <FONT color="green">227</FONT> getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {<a name="line.227"></a>
231 <FONT color="green">228</FONT> /** {@inheritDoc} */<a name="line.228"></a>
232 <FONT color="green">229</FONT> @Override<a name="line.229"></a>
233 <FONT color="green">230</FONT> public void visit(final int row, final int column,<a name="line.230"></a>
234 <FONT color="green">231</FONT> final double value) {<a name="line.231"></a>
235 <FONT color="green">232</FONT> data[row][column] = value / singularValues[row];<a name="line.232"></a>
236 <FONT color="green">233</FONT> }<a name="line.233"></a>
237 <FONT color="green">234</FONT> }, 0, dimension - 1, 0, p - 1);<a name="line.234"></a>
238 <FONT color="green">235</FONT> <a name="line.235"></a>
239 <FONT color="green">236</FONT> RealMatrix jv = new Array2DRowRealMatrix(data, false);<a name="line.236"></a>
240 <FONT color="green">237</FONT> return jv.transpose().multiply(jv);<a name="line.237"></a>
241 <FONT color="green">238</FONT> <a name="line.238"></a>
242 <FONT color="green">239</FONT> }<a name="line.239"></a>
243 <FONT color="green">240</FONT> <a name="line.240"></a>
244 <FONT color="green">241</FONT> /** {@inheritDoc} */<a name="line.241"></a>
245 <FONT color="green">242</FONT> public double getNorm() throws InvalidMatrixException {<a name="line.242"></a>
246 <FONT color="green">243</FONT> return singularValues[0];<a name="line.243"></a>
247 <FONT color="green">244</FONT> }<a name="line.244"></a>
248 <FONT color="green">245</FONT> <a name="line.245"></a>
249 <FONT color="green">246</FONT> /** {@inheritDoc} */<a name="line.246"></a>
250 <FONT color="green">247</FONT> public double getConditionNumber() throws InvalidMatrixException {<a name="line.247"></a>
251 <FONT color="green">248</FONT> return singularValues[0] / singularValues[singularValues.length - 1];<a name="line.248"></a>
252 <FONT color="green">249</FONT> }<a name="line.249"></a>
253 <FONT color="green">250</FONT> <a name="line.250"></a>
254 <FONT color="green">251</FONT> /** {@inheritDoc} */<a name="line.251"></a>
255 <FONT color="green">252</FONT> public int getRank() throws IllegalStateException {<a name="line.252"></a>
256 <FONT color="green">253</FONT> <a name="line.253"></a>
257 <FONT color="green">254</FONT> final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);<a name="line.254"></a>
258 <FONT color="green">255</FONT> <a name="line.255"></a>
259 <FONT color="green">256</FONT> for (int i = singularValues.length - 1; i &gt;= 0; --i) {<a name="line.256"></a>
260 <FONT color="green">257</FONT> if (singularValues[i] &gt; threshold) {<a name="line.257"></a>
261 <FONT color="green">258</FONT> return i + 1;<a name="line.258"></a>
262 <FONT color="green">259</FONT> }<a name="line.259"></a>
263 <FONT color="green">260</FONT> }<a name="line.260"></a>
264 <FONT color="green">261</FONT> return 0;<a name="line.261"></a>
265 <FONT color="green">262</FONT> <a name="line.262"></a>
266 <FONT color="green">263</FONT> }<a name="line.263"></a>
267 <FONT color="green">264</FONT> <a name="line.264"></a>
268 <FONT color="green">265</FONT> /** {@inheritDoc} */<a name="line.265"></a>
269 <FONT color="green">266</FONT> public DecompositionSolver getSolver() {<a name="line.266"></a>
270 <FONT color="green">267</FONT> return new Solver(singularValues, getUT(), getV(), getRank() == Math<a name="line.267"></a>
271 <FONT color="green">268</FONT> .max(m, n));<a name="line.268"></a>
272 <FONT color="green">269</FONT> }<a name="line.269"></a>
273 <FONT color="green">270</FONT> <a name="line.270"></a>
274 <FONT color="green">271</FONT> /** Specialized solver. */<a name="line.271"></a>
275 <FONT color="green">272</FONT> private static class Solver implements DecompositionSolver {<a name="line.272"></a>
276 <FONT color="green">273</FONT> <a name="line.273"></a>
277 <FONT color="green">274</FONT> /** Pseudo-inverse of the initial matrix. */<a name="line.274"></a>
278 <FONT color="green">275</FONT> private final RealMatrix pseudoInverse;<a name="line.275"></a>
279 <FONT color="green">276</FONT> <a name="line.276"></a>
280 <FONT color="green">277</FONT> /** Singularity indicator. */<a name="line.277"></a>
281 <FONT color="green">278</FONT> private boolean nonSingular;<a name="line.278"></a>
282 <FONT color="green">279</FONT> <a name="line.279"></a>
283 <FONT color="green">280</FONT> /**<a name="line.280"></a>
284 <FONT color="green">281</FONT> * Build a solver from decomposed matrix.<a name="line.281"></a>
285 <FONT color="green">282</FONT> * @param singularValues<a name="line.282"></a>
286 <FONT color="green">283</FONT> * singularValues<a name="line.283"></a>
287 <FONT color="green">284</FONT> * @param uT<a name="line.284"></a>
288 <FONT color="green">285</FONT> * U&lt;sup&gt;T&lt;/sup&gt; matrix of the decomposition<a name="line.285"></a>
289 <FONT color="green">286</FONT> * @param v<a name="line.286"></a>
290 <FONT color="green">287</FONT> * V matrix of the decomposition<a name="line.287"></a>
291 <FONT color="green">288</FONT> * @param nonSingular<a name="line.288"></a>
292 <FONT color="green">289</FONT> * singularity indicator<a name="line.289"></a>
293 <FONT color="green">290</FONT> */<a name="line.290"></a>
294 <FONT color="green">291</FONT> private Solver(final double[] singularValues, final RealMatrix uT,<a name="line.291"></a>
295 <FONT color="green">292</FONT> final RealMatrix v, final boolean nonSingular) {<a name="line.292"></a>
296 <FONT color="green">293</FONT> double[][] suT = uT.getData();<a name="line.293"></a>
297 <FONT color="green">294</FONT> for (int i = 0; i &lt; singularValues.length; ++i) {<a name="line.294"></a>
298 <FONT color="green">295</FONT> final double a;<a name="line.295"></a>
299 <FONT color="green">296</FONT> if (singularValues[i]&gt;0) {<a name="line.296"></a>
300 <FONT color="green">297</FONT> a=1.0 / singularValues[i];<a name="line.297"></a>
301 <FONT color="green">298</FONT> } else {<a name="line.298"></a>
302 <FONT color="green">299</FONT> a=0.0;<a name="line.299"></a>
303 <FONT color="green">300</FONT> }<a name="line.300"></a>
304 <FONT color="green">301</FONT> final double[] suTi = suT[i];<a name="line.301"></a>
305 <FONT color="green">302</FONT> for (int j = 0; j &lt; suTi.length; ++j) {<a name="line.302"></a>
306 <FONT color="green">303</FONT> suTi[j] *= a;<a name="line.303"></a>
307 <FONT color="green">304</FONT> }<a name="line.304"></a>
308 <FONT color="green">305</FONT> }<a name="line.305"></a>
309 <FONT color="green">306</FONT> pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));<a name="line.306"></a>
310 <FONT color="green">307</FONT> this.nonSingular = nonSingular;<a name="line.307"></a>
311 <FONT color="green">308</FONT> }<a name="line.308"></a>
312 <FONT color="green">309</FONT> <a name="line.309"></a>
313 <FONT color="green">310</FONT> /**<a name="line.310"></a>
314 <FONT color="green">311</FONT> * Solve the linear equation A &amp;times; X = B in least square sense.<a name="line.311"></a>
315 <FONT color="green">312</FONT> * &lt;p&gt;<a name="line.312"></a>
316 <FONT color="green">313</FONT> * The m&amp;times;n matrix A may not be square, the solution X is such that<a name="line.313"></a>
317 <FONT color="green">314</FONT> * ||A &amp;times; X - B|| is minimal.<a name="line.314"></a>
318 <FONT color="green">315</FONT> * &lt;/p&gt;<a name="line.315"></a>
319 <FONT color="green">316</FONT> * @param b<a name="line.316"></a>
320 <FONT color="green">317</FONT> * right-hand side of the equation A &amp;times; X = B<a name="line.317"></a>
321 <FONT color="green">318</FONT> * @return a vector X that minimizes the two norm of A &amp;times; X - B<a name="line.318"></a>
322 <FONT color="green">319</FONT> * @exception IllegalArgumentException<a name="line.319"></a>
323 <FONT color="green">320</FONT> * if matrices dimensions don't match<a name="line.320"></a>
324 <FONT color="green">321</FONT> */<a name="line.321"></a>
325 <FONT color="green">322</FONT> public double[] solve(final double[] b) throws IllegalArgumentException {<a name="line.322"></a>
326 <FONT color="green">323</FONT> return pseudoInverse.operate(b);<a name="line.323"></a>
327 <FONT color="green">324</FONT> }<a name="line.324"></a>
328 <FONT color="green">325</FONT> <a name="line.325"></a>
329 <FONT color="green">326</FONT> /**<a name="line.326"></a>
330 <FONT color="green">327</FONT> * Solve the linear equation A &amp;times; X = B in least square sense.<a name="line.327"></a>
331 <FONT color="green">328</FONT> * &lt;p&gt;<a name="line.328"></a>
332 <FONT color="green">329</FONT> * The m&amp;times;n matrix A may not be square, the solution X is such that<a name="line.329"></a>
333 <FONT color="green">330</FONT> * ||A &amp;times; X - B|| is minimal.<a name="line.330"></a>
334 <FONT color="green">331</FONT> * &lt;/p&gt;<a name="line.331"></a>
335 <FONT color="green">332</FONT> * @param b<a name="line.332"></a>
336 <FONT color="green">333</FONT> * right-hand side of the equation A &amp;times; X = B<a name="line.333"></a>
337 <FONT color="green">334</FONT> * @return a vector X that minimizes the two norm of A &amp;times; X - B<a name="line.334"></a>
338 <FONT color="green">335</FONT> * @exception IllegalArgumentException<a name="line.335"></a>
339 <FONT color="green">336</FONT> * if matrices dimensions don't match<a name="line.336"></a>
340 <FONT color="green">337</FONT> */<a name="line.337"></a>
341 <FONT color="green">338</FONT> public RealVector solve(final RealVector b)<a name="line.338"></a>
342 <FONT color="green">339</FONT> throws IllegalArgumentException {<a name="line.339"></a>
343 <FONT color="green">340</FONT> return pseudoInverse.operate(b);<a name="line.340"></a>
344 <FONT color="green">341</FONT> }<a name="line.341"></a>
345 <FONT color="green">342</FONT> <a name="line.342"></a>
346 <FONT color="green">343</FONT> /**<a name="line.343"></a>
347 <FONT color="green">344</FONT> * Solve the linear equation A &amp;times; X = B in least square sense.<a name="line.344"></a>
348 <FONT color="green">345</FONT> * &lt;p&gt;<a name="line.345"></a>
349 <FONT color="green">346</FONT> * The m&amp;times;n matrix A may not be square, the solution X is such that<a name="line.346"></a>
350 <FONT color="green">347</FONT> * ||A &amp;times; X - B|| is minimal.<a name="line.347"></a>
351 <FONT color="green">348</FONT> * &lt;/p&gt;<a name="line.348"></a>
352 <FONT color="green">349</FONT> * @param b<a name="line.349"></a>
353 <FONT color="green">350</FONT> * right-hand side of the equation A &amp;times; X = B<a name="line.350"></a>
354 <FONT color="green">351</FONT> * @return a matrix X that minimizes the two norm of A &amp;times; X - B<a name="line.351"></a>
355 <FONT color="green">352</FONT> * @exception IllegalArgumentException<a name="line.352"></a>
356 <FONT color="green">353</FONT> * if matrices dimensions don't match<a name="line.353"></a>
357 <FONT color="green">354</FONT> */<a name="line.354"></a>
358 <FONT color="green">355</FONT> public RealMatrix solve(final RealMatrix b)<a name="line.355"></a>
359 <FONT color="green">356</FONT> throws IllegalArgumentException {<a name="line.356"></a>
360 <FONT color="green">357</FONT> return pseudoInverse.multiply(b);<a name="line.357"></a>
361 <FONT color="green">358</FONT> }<a name="line.358"></a>
362 <FONT color="green">359</FONT> <a name="line.359"></a>
363 <FONT color="green">360</FONT> /**<a name="line.360"></a>
364 <FONT color="green">361</FONT> * Check if the decomposed matrix is non-singular.<a name="line.361"></a>
365 <FONT color="green">362</FONT> * @return true if the decomposed matrix is non-singular<a name="line.362"></a>
366 <FONT color="green">363</FONT> */<a name="line.363"></a>
367 <FONT color="green">364</FONT> public boolean isNonSingular() {<a name="line.364"></a>
368 <FONT color="green">365</FONT> return nonSingular;<a name="line.365"></a>
369 <FONT color="green">366</FONT> }<a name="line.366"></a>
370 <FONT color="green">367</FONT> <a name="line.367"></a>
371 <FONT color="green">368</FONT> /**<a name="line.368"></a>
372 <FONT color="green">369</FONT> * Get the pseudo-inverse of the decomposed matrix.<a name="line.369"></a>
373 <FONT color="green">370</FONT> * @return inverse matrix<a name="line.370"></a>
374 <FONT color="green">371</FONT> */<a name="line.371"></a>
375 <FONT color="green">372</FONT> public RealMatrix getInverse() {<a name="line.372"></a>
376 <FONT color="green">373</FONT> return pseudoInverse;<a name="line.373"></a>
377 <FONT color="green">374</FONT> }<a name="line.374"></a>
378 <FONT color="green">375</FONT> <a name="line.375"></a>
379 <FONT color="green">376</FONT> }<a name="line.376"></a>
380 <FONT color="green">377</FONT> <a name="line.377"></a>
381 <FONT color="green">378</FONT> }<a name="line.378"></a>
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442 </PRE>
443 </BODY>
444 </HTML>