Free cookie consent management tool by TermsFeed Policy Generator

source: tags/3.3.2/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/svd.cs @ 13398

Last change on this file since 13398 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)

File size: 14.9 KB
Line 
1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class svd
26    {
27        /*************************************************************************
28        Singular value decomposition of a rectangular matrix.
29
30        The algorithm calculates the singular value decomposition of a matrix of
31        size MxN: A = U * S * V^T
32
33        The algorithm finds the singular values and, optionally, matrices U and V^T.
34        The algorithm can find both first min(M,N) columns of matrix U and rows of
35        matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
36        and NxN respectively).
37
38        Take into account that the subroutine does not return matrix V but V^T.
39
40        Input parameters:
41            A           -   matrix to be decomposed.
42                            Array whose indexes range within [0..M-1, 0..N-1].
43            M           -   number of rows in matrix A.
44            N           -   number of columns in matrix A.
45            UNeeded     -   0, 1 or 2. See the description of the parameter U.
46            VTNeeded    -   0, 1 or 2. See the description of the parameter VT.
47            AdditionalMemory -
48                            If the parameter:
49                             * equals 0, the algorithm doesn’t use additional
50                               memory (lower requirements, lower performance).
51                             * equals 1, the algorithm uses additional
52                               memory of size min(M,N)*min(M,N) of real numbers.
53                               It often speeds up the algorithm.
54                             * equals 2, the algorithm uses additional
55                               memory of size M*min(M,N) of real numbers.
56                               It allows to get a maximum performance.
57                            The recommended value of the parameter is 2.
58
59        Output parameters:
60            W           -   contains singular values in descending order.
61            U           -   if UNeeded=0, U isn't changed, the left singular vectors
62                            are not calculated.
63                            if Uneeded=1, U contains left singular vectors (first
64                            min(M,N) columns of matrix U). Array whose indexes range
65                            within [0..M-1, 0..Min(M,N)-1].
66                            if UNeeded=2, U contains matrix U wholly. Array whose
67                            indexes range within [0..M-1, 0..M-1].
68            VT          -   if VTNeeded=0, VT isn’t changed, the right singular vectors
69                            are not calculated.
70                            if VTNeeded=1, VT contains right singular vectors (first
71                            min(M,N) rows of matrix V^T). Array whose indexes range
72                            within [0..min(M,N)-1, 0..N-1].
73                            if VTNeeded=2, VT contains matrix V^T wholly. Array whose
74                            indexes range within [0..N-1, 0..N-1].
75
76          -- ALGLIB --
77             Copyright 2005 by Bochkanov Sergey
78        *************************************************************************/
79        public static bool rmatrixsvd(double[,] a,
80            int m,
81            int n,
82            int uneeded,
83            int vtneeded,
84            int additionalmemory,
85            ref double[] w,
86            ref double[,] u,
87            ref double[,] vt)
88        {
89            bool result = new bool();
90            double[] tauq = new double[0];
91            double[] taup = new double[0];
92            double[] tau = new double[0];
93            double[] e = new double[0];
94            double[] work = new double[0];
95            double[,] t2 = new double[0,0];
96            bool isupper = new bool();
97            int minmn = 0;
98            int ncu = 0;
99            int nrvt = 0;
100            int nru = 0;
101            int ncvt = 0;
102            int i = 0;
103            int j = 0;
104
105            a = (double[,])a.Clone();
106
107            result = true;
108            if( m==0 | n==0 )
109            {
110                return result;
111            }
112            System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
113            System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
114            System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
115           
116            //
117            // initialize
118            //
119            minmn = Math.Min(m, n);
120            w = new double[minmn+1];
121            ncu = 0;
122            nru = 0;
123            if( uneeded==1 )
124            {
125                nru = m;
126                ncu = minmn;
127                u = new double[nru-1+1, ncu-1+1];
128            }
129            if( uneeded==2 )
130            {
131                nru = m;
132                ncu = m;
133                u = new double[nru-1+1, ncu-1+1];
134            }
135            nrvt = 0;
136            ncvt = 0;
137            if( vtneeded==1 )
138            {
139                nrvt = minmn;
140                ncvt = n;
141                vt = new double[nrvt-1+1, ncvt-1+1];
142            }
143            if( vtneeded==2 )
144            {
145                nrvt = n;
146                ncvt = n;
147                vt = new double[nrvt-1+1, ncvt-1+1];
148            }
149           
150            //
151            // M much larger than N
152            // Use bidiagonal reduction with QR-decomposition
153            //
154            if( (double)(m)>(double)(1.6*n) )
155            {
156                if( uneeded==0 )
157                {
158                   
159                    //
160                    // No left singular vectors to be computed
161                    //
162                    ortfac.rmatrixqr(ref a, m, n, ref tau);
163                    for(i=0; i<=n-1; i++)
164                    {
165                        for(j=0; j<=i-1; j++)
166                        {
167                            a[i,j] = 0;
168                        }
169                    }
170                    ortfac.rmatrixbd(ref a, n, n, ref tauq, ref taup);
171                    ortfac.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
172                    ortfac.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
173                    result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
174                    return result;
175                }
176                else
177                {
178                   
179                    //
180                    // Left singular vectors (may be full matrix U) to be computed
181                    //
182                    ortfac.rmatrixqr(ref a, m, n, ref tau);
183                    ortfac.rmatrixqrunpackq(ref a, m, n, ref tau, ncu, ref u);
184                    for(i=0; i<=n-1; i++)
185                    {
186                        for(j=0; j<=i-1; j++)
187                        {
188                            a[i,j] = 0;
189                        }
190                    }
191                    ortfac.rmatrixbd(ref a, n, n, ref tauq, ref taup);
192                    ortfac.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
193                    ortfac.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
194                    if( additionalmemory<1 )
195                    {
196                       
197                        //
198                        // No additional memory can be used
199                        //
200                        ortfac.rmatrixbdmultiplybyq(ref a, n, n, ref tauq, ref u, m, n, true, false);
201                        result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
202                    }
203                    else
204                    {
205                       
206                        //
207                        // Large U. Transforming intermediate matrix T2
208                        //
209                        work = new double[Math.Max(m, n)+1];
210                        ortfac.rmatrixbdunpackq(ref a, n, n, ref tauq, n, ref t2);
211                        blas.copymatrix(ref u, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
212                        blas.inplacetranspose(ref t2, 0, n-1, 0, n-1, ref work);
213                        result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
214                        blas.matrixmatrixmultiply(ref a, 0, m-1, 0, n-1, false, ref t2, 0, n-1, 0, n-1, true, 1.0, ref u, 0, m-1, 0, n-1, 0.0, ref work);
215                    }
216                    return result;
217                }
218            }
219           
220            //
221            // N much larger than M
222            // Use bidiagonal reduction with LQ-decomposition
223            //
224            if( (double)(n)>(double)(1.6*m) )
225            {
226                if( vtneeded==0 )
227                {
228                   
229                    //
230                    // No right singular vectors to be computed
231                    //
232                    ortfac.rmatrixlq(ref a, m, n, ref tau);
233                    for(i=0; i<=m-1; i++)
234                    {
235                        for(j=i+1; j<=m-1; j++)
236                        {
237                            a[i,j] = 0;
238                        }
239                    }
240                    ortfac.rmatrixbd(ref a, m, m, ref tauq, ref taup);
241                    ortfac.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
242                    ortfac.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
243                    work = new double[m+1];
244                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
245                    result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
246                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
247                    return result;
248                }
249                else
250                {
251                   
252                    //
253                    // Right singular vectors (may be full matrix VT) to be computed
254                    //
255                    ortfac.rmatrixlq(ref a, m, n, ref tau);
256                    ortfac.rmatrixlqunpackq(ref a, m, n, ref tau, nrvt, ref vt);
257                    for(i=0; i<=m-1; i++)
258                    {
259                        for(j=i+1; j<=m-1; j++)
260                        {
261                            a[i,j] = 0;
262                        }
263                    }
264                    ortfac.rmatrixbd(ref a, m, m, ref tauq, ref taup);
265                    ortfac.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
266                    ortfac.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
267                    work = new double[Math.Max(m, n)+1];
268                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
269                    if( additionalmemory<1 )
270                    {
271                       
272                        //
273                        // No additional memory available
274                        //
275                        ortfac.rmatrixbdmultiplybyp(ref a, m, m, ref taup, ref vt, m, n, false, true);
276                        result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
277                    }
278                    else
279                    {
280                       
281                        //
282                        // Large VT. Transforming intermediate matrix T2
283                        //
284                        ortfac.rmatrixbdunpackpt(ref a, m, m, ref taup, m, ref t2);
285                        result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
286                        blas.copymatrix(ref vt, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
287                        blas.matrixmatrixmultiply(ref t2, 0, m-1, 0, m-1, false, ref a, 0, m-1, 0, n-1, false, 1.0, ref vt, 0, m-1, 0, n-1, 0.0, ref work);
288                    }
289                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
290                    return result;
291                }
292            }
293           
294            //
295            // M<=N
296            // We can use inplace transposition of U to get rid of columnwise operations
297            //
298            if( m<=n )
299            {
300                ortfac.rmatrixbd(ref a, m, n, ref tauq, ref taup);
301                ortfac.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
302                ortfac.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
303                ortfac.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
304                work = new double[m+1];
305                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
306                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
307                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
308                return result;
309            }
310           
311            //
312            // Simple bidiagonal reduction
313            //
314            ortfac.rmatrixbd(ref a, m, n, ref tauq, ref taup);
315            ortfac.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
316            ortfac.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
317            ortfac.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
318            if( additionalmemory<2 | uneeded==0 )
319            {
320               
321                //
322                // We cant use additional memory or there is no need in such operations
323                //
324                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
325            }
326            else
327            {
328               
329                //
330                // We can use additional memory
331                //
332                t2 = new double[minmn-1+1, m-1+1];
333                blas.copyandtranspose(ref u, 0, m-1, 0, minmn-1, ref t2, 0, minmn-1, 0, m-1);
334                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
335                blas.copyandtranspose(ref t2, 0, minmn-1, 0, m-1, ref u, 0, m-1, 0, minmn-1);
336            }
337            return result;
338        }
339    }
340}
Note: See TracBrowser for help on using the repository browser.