Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.LinearRegression/3.2/alglib/qr.cs @ 2211

Last change on this file since 2211 was 2154, checked in by gkronber, 15 years ago

Added linear regression plugin. #697

File size: 13.7 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10Redistribution and use in source and binary forms, with or without
11modification, are permitted provided that the following conditions are
12met:
13
14- Redistributions of source code must retain the above copyright
15  notice, this list of conditions and the following disclaimer.
16
17- Redistributions in binary form must reproduce the above copyright
18  notice, this list of conditions and the following disclaimer listed
19  in this license in the documentation and/or other materials
20  provided with the distribution.
21
22- Neither the name of the copyright holders nor the names of its
23  contributors may be used to endorse or promote products derived from
24  this software without specific prior written permission.
25
26THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37*************************************************************************/
38
39using System;
40
41class qr
42{
43    /*************************************************************************
44    QR decomposition of a rectangular matrix of size MxN
45
46    Input parameters:
47        A   -   matrix A whose indexes range within [0..M-1, 0..N-1].
48        M   -   number of rows in matrix A.
49        N   -   number of columns in matrix A.
50
51    Output parameters:
52        A   -   matrices Q and R in compact form (see below).
53        Tau -   array of scalar factors which are used to form
54                matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)].
55
56    Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
57    MxM, R - upper triangular (or upper trapezoid) matrix of size M x N.
58
59    The elements of matrix R are located on and above the main diagonal of
60    matrix A. The elements which are located in Tau array and below the main
61    diagonal of matrix A are used to form matrix Q as follows:
62
63    Matrix Q is represented as a product of elementary reflections
64
65    Q = H(0)*H(2)*...*H(k-1),
66
67    where k = min(m,n), and each H(i) is in the form
68
69    H(i) = 1 - tau * v * (v^T)
70
71    where tau is a scalar stored in Tau[I]; v - real vector,
72    so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i).
73
74      -- LAPACK routine (version 3.0) --
75         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
76         Courant Institute, Argonne National Lab, and Rice University
77         February 29, 1992.
78         Translation from FORTRAN to pseudocode (AlgoPascal)
79         by Sergey Bochkanov, ALGLIB project, 2005-2007.
80    *************************************************************************/
81    public static void rmatrixqr(ref double[,] a,
82        int m,
83        int n,
84        ref double[] tau)
85    {
86        double[] work = new double[0];
87        double[] t = new double[0];
88        int i = 0;
89        int k = 0;
90        int minmn = 0;
91        double tmp = 0;
92        int i_ = 0;
93        int i1_ = 0;
94
95        if( m<=0 | n<=0 )
96        {
97            return;
98        }
99        minmn = Math.Min(m, n);
100        work = new double[n-1+1];
101        t = new double[m+1];
102        tau = new double[minmn-1+1];
103       
104        //
105        // Test the input arguments
106        //
107        k = minmn;
108        for(i=0; i<=k-1; i++)
109        {
110           
111            //
112            // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
113            //
114            i1_ = (i) - (1);
115            for(i_=1; i_<=m-i;i_++)
116            {
117                t[i_] = a[i_+i1_,i];
118            }
119            reflections.generatereflection(ref t, m-i, ref tmp);
120            tau[i] = tmp;
121            i1_ = (1) - (i);
122            for(i_=i; i_<=m-1;i_++)
123            {
124                a[i_,i] = t[i_+i1_];
125            }
126            t[1] = 1;
127            if( i<n )
128            {
129               
130                //
131                // Apply H(i) to A(i:m-1,i+1:n-1) from the left
132                //
133                reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m-1, i+1, n-1, ref work);
134            }
135        }
136    }
137
138
139    /*************************************************************************
140    Partial unpacking of matrix Q from the QR decomposition of a matrix A
141
142    Input parameters:
143        A       -   matrices Q and R in compact form.
144                    Output of RMatrixQR subroutine.
145        M       -   number of rows in given matrix A. M>=0.
146        N       -   number of columns in given matrix A. N>=0.
147        Tau     -   scalar factors which are used to form Q.
148                    Output of the RMatrixQR subroutine.
149        QColumns -  required number of columns of matrix Q. M>=QColumns>=0.
150
151    Output parameters:
152        Q       -   first QColumns columns of matrix Q.
153                    Array whose indexes range within [0..M-1, 0..QColumns-1].
154                    If QColumns=0, the array remains unchanged.
155
156      -- ALGLIB --
157         Copyright 2005 by Bochkanov Sergey
158    *************************************************************************/
159    public static void rmatrixqrunpackq(ref double[,] a,
160        int m,
161        int n,
162        ref double[] tau,
163        int qcolumns,
164        ref double[,] q)
165    {
166        int i = 0;
167        int j = 0;
168        int k = 0;
169        int minmn = 0;
170        double[] v = new double[0];
171        double[] work = new double[0];
172        int i_ = 0;
173        int i1_ = 0;
174
175        System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
176        if( m<=0 | n<=0 | qcolumns<=0 )
177        {
178            return;
179        }
180       
181        //
182        // init
183        //
184        minmn = Math.Min(m, n);
185        k = Math.Min(minmn, qcolumns);
186        q = new double[m-1+1, qcolumns-1+1];
187        v = new double[m+1];
188        work = new double[qcolumns-1+1];
189        for(i=0; i<=m-1; i++)
190        {
191            for(j=0; j<=qcolumns-1; j++)
192            {
193                if( i==j )
194                {
195                    q[i,j] = 1;
196                }
197                else
198                {
199                    q[i,j] = 0;
200                }
201            }
202        }
203       
204        //
205        // unpack Q
206        //
207        for(i=k-1; i>=0; i--)
208        {
209           
210            //
211            // Apply H(i)
212            //
213            i1_ = (i) - (1);
214            for(i_=1; i_<=m-i;i_++)
215            {
216                v[i_] = a[i_+i1_,i];
217            }
218            v[1] = 1;
219            reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m-1, 0, qcolumns-1, ref work);
220        }
221    }
222
223
224    /*************************************************************************
225    Unpacking of matrix R from the QR decomposition of a matrix A
226
227    Input parameters:
228        A       -   matrices Q and R in compact form.
229                    Output of RMatrixQR subroutine.
230        M       -   number of rows in given matrix A. M>=0.
231        N       -   number of columns in given matrix A. N>=0.
232
233    Output parameters:
234        R       -   matrix R, array[0..M-1, 0..N-1].
235
236      -- ALGLIB --
237         Copyright 2005 by Bochkanov Sergey
238    *************************************************************************/
239    public static void rmatrixqrunpackr(ref double[,] a,
240        int m,
241        int n,
242        ref double[,] r)
243    {
244        int i = 0;
245        int k = 0;
246        int i_ = 0;
247
248        if( m<=0 | n<=0 )
249        {
250            return;
251        }
252        k = Math.Min(m, n);
253        r = new double[m-1+1, n-1+1];
254        for(i=0; i<=n-1; i++)
255        {
256            r[0,i] = 0;
257        }
258        for(i=1; i<=m-1; i++)
259        {
260            for(i_=0; i_<=n-1;i_++)
261            {
262                r[i,i_] = r[0,i_];
263            }
264        }
265        for(i=0; i<=k-1; i++)
266        {
267            for(i_=i; i_<=n-1;i_++)
268            {
269                r[i,i_] = a[i,i_];
270            }
271        }
272    }
273
274
275    /*************************************************************************
276    Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
277    *************************************************************************/
278    public static void qrdecomposition(ref double[,] a,
279        int m,
280        int n,
281        ref double[] tau)
282    {
283        double[] work = new double[0];
284        double[] t = new double[0];
285        int i = 0;
286        int k = 0;
287        int mmip1 = 0;
288        int minmn = 0;
289        double tmp = 0;
290        int i_ = 0;
291        int i1_ = 0;
292
293        minmn = Math.Min(m, n);
294        work = new double[n+1];
295        t = new double[m+1];
296        tau = new double[minmn+1];
297       
298        //
299        // Test the input arguments
300        //
301        k = Math.Min(m, n);
302        for(i=1; i<=k; i++)
303        {
304           
305            //
306            // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
307            //
308            mmip1 = m-i+1;
309            i1_ = (i) - (1);
310            for(i_=1; i_<=mmip1;i_++)
311            {
312                t[i_] = a[i_+i1_,i];
313            }
314            reflections.generatereflection(ref t, mmip1, ref tmp);
315            tau[i] = tmp;
316            i1_ = (1) - (i);
317            for(i_=i; i_<=m;i_++)
318            {
319                a[i_,i] = t[i_+i1_];
320            }
321            t[1] = 1;
322            if( i<n )
323            {
324               
325                //
326                // Apply H(i) to A(i:m,i+1:n) from the left
327                //
328                reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m, i+1, n, ref work);
329            }
330        }
331    }
332
333
334    /*************************************************************************
335    Obsolete 1-based subroutine. See RMatrixQRUnpackQ for 0-based replacement.
336    *************************************************************************/
337    public static void unpackqfromqr(ref double[,] a,
338        int m,
339        int n,
340        ref double[] tau,
341        int qcolumns,
342        ref double[,] q)
343    {
344        int i = 0;
345        int j = 0;
346        int k = 0;
347        int minmn = 0;
348        double[] v = new double[0];
349        double[] work = new double[0];
350        int vm = 0;
351        int i_ = 0;
352        int i1_ = 0;
353
354        System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
355        if( m==0 | n==0 | qcolumns==0 )
356        {
357            return;
358        }
359       
360        //
361        // init
362        //
363        minmn = Math.Min(m, n);
364        k = Math.Min(minmn, qcolumns);
365        q = new double[m+1, qcolumns+1];
366        v = new double[m+1];
367        work = new double[qcolumns+1];
368        for(i=1; i<=m; i++)
369        {
370            for(j=1; j<=qcolumns; j++)
371            {
372                if( i==j )
373                {
374                    q[i,j] = 1;
375                }
376                else
377                {
378                    q[i,j] = 0;
379                }
380            }
381        }
382       
383        //
384        // unpack Q
385        //
386        for(i=k; i>=1; i--)
387        {
388           
389            //
390            // Apply H(i)
391            //
392            vm = m-i+1;
393            i1_ = (i) - (1);
394            for(i_=1; i_<=vm;i_++)
395            {
396                v[i_] = a[i_+i1_,i];
397            }
398            v[1] = 1;
399            reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m, 1, qcolumns, ref work);
400        }
401    }
402
403
404    /*************************************************************************
405    Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
406    *************************************************************************/
407    public static void qrdecompositionunpacked(double[,] a,
408        int m,
409        int n,
410        ref double[,] q,
411        ref double[,] r)
412    {
413        int i = 0;
414        int k = 0;
415        double[] tau = new double[0];
416        double[] work = new double[0];
417        double[] v = new double[0];
418        int i_ = 0;
419
420        a = (double[,])a.Clone();
421
422        k = Math.Min(m, n);
423        if( n<=0 )
424        {
425            return;
426        }
427        work = new double[m+1];
428        v = new double[m+1];
429        q = new double[m+1, m+1];
430        r = new double[m+1, n+1];
431       
432        //
433        // QRDecomposition
434        //
435        qrdecomposition(ref a, m, n, ref tau);
436       
437        //
438        // R
439        //
440        for(i=1; i<=n; i++)
441        {
442            r[1,i] = 0;
443        }
444        for(i=2; i<=m; i++)
445        {
446            for(i_=1; i_<=n;i_++)
447            {
448                r[i,i_] = r[1,i_];
449            }
450        }
451        for(i=1; i<=k; i++)
452        {
453            for(i_=i; i_<=n;i_++)
454            {
455                r[i,i_] = a[i,i_];
456            }
457        }
458       
459        //
460        // Q
461        //
462        unpackqfromqr(ref a, m, n, ref tau, m, ref q);
463    }
464}
Note: See TracBrowser for help on using the repository browser.