Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/inverseupdate.cs @ 2643

Last change on this file since 2643 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 15.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 inverseupdate
26    {
27        /*************************************************************************
28        Inverse matrix update by the Sherman-Morrison formula
29
30        The algorithm updates matrix A^-1 when adding a number to an element
31        of matrix A.
32
33        Input parameters:
34            InvA    -   inverse of matrix A.
35                        Array whose indexes range within [0..N-1, 0..N-1].
36            N       -   size of matrix A.
37            UpdRow  -   row where the element to be updated is stored.
38            UpdColumn - column where the element to be updated is stored.
39            UpdVal  -   a number to be added to the element.
40
41
42        Output parameters:
43            InvA    -   inverse of modified matrix A.
44
45          -- ALGLIB --
46             Copyright 2005 by Bochkanov Sergey
47        *************************************************************************/
48        public static void rmatrixinvupdatesimple(ref double[,] inva,
49            int n,
50            int updrow,
51            int updcolumn,
52            double updval)
53        {
54            double[] t1 = new double[0];
55            double[] t2 = new double[0];
56            int i = 0;
57            int j = 0;
58            double lambda = 0;
59            double vt = 0;
60            int i_ = 0;
61
62            System.Diagnostics.Debug.Assert(updrow>=0 & updrow<n, "RMatrixInvUpdateSimple: incorrect UpdRow!");
63            System.Diagnostics.Debug.Assert(updcolumn>=0 & updcolumn<n, "RMatrixInvUpdateSimple: incorrect UpdColumn!");
64            t1 = new double[n-1+1];
65            t2 = new double[n-1+1];
66           
67            //
68            // T1 = InvA * U
69            //
70            for(i_=0; i_<=n-1;i_++)
71            {
72                t1[i_] = inva[i_,updrow];
73            }
74           
75            //
76            // T2 = v*InvA
77            //
78            for(i_=0; i_<=n-1;i_++)
79            {
80                t2[i_] = inva[updcolumn,i_];
81            }
82           
83            //
84            // Lambda = v * InvA * U
85            //
86            lambda = updval*inva[updcolumn,updrow];
87           
88            //
89            // InvA = InvA - correction
90            //
91            for(i=0; i<=n-1; i++)
92            {
93                vt = updval*t1[i];
94                vt = vt/(1+lambda);
95                for(i_=0; i_<=n-1;i_++)
96                {
97                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
98                }
99            }
100        }
101
102
103        /*************************************************************************
104        Inverse matrix update by the Sherman-Morrison formula
105
106        The algorithm updates matrix A^-1 when adding a vector to a row
107        of matrix A.
108
109        Input parameters:
110            InvA    -   inverse of matrix A.
111                        Array whose indexes range within [0..N-1, 0..N-1].
112            N       -   size of matrix A.
113            UpdRow  -   the row of A whose vector V was added.
114                        0 <= Row <= N-1
115            V       -   the vector to be added to a row.
116                        Array whose index ranges within [0..N-1].
117
118        Output parameters:
119            InvA    -   inverse of modified matrix A.
120
121          -- ALGLIB --
122             Copyright 2005 by Bochkanov Sergey
123        *************************************************************************/
124        public static void rmatrixinvupdaterow(ref double[,] inva,
125            int n,
126            int updrow,
127            ref double[] v)
128        {
129            double[] t1 = new double[0];
130            double[] t2 = new double[0];
131            int i = 0;
132            int j = 0;
133            double lambda = 0;
134            double vt = 0;
135            int i_ = 0;
136
137            t1 = new double[n-1+1];
138            t2 = new double[n-1+1];
139           
140            //
141            // T1 = InvA * U
142            //
143            for(i_=0; i_<=n-1;i_++)
144            {
145                t1[i_] = inva[i_,updrow];
146            }
147           
148            //
149            // T2 = v*InvA
150            // Lambda = v * InvA * U
151            //
152            for(j=0; j<=n-1; j++)
153            {
154                vt = 0.0;
155                for(i_=0; i_<=n-1;i_++)
156                {
157                    vt += v[i_]*inva[i_,j];
158                }
159                t2[j] = vt;
160            }
161            lambda = t2[updrow];
162           
163            //
164            // InvA = InvA - correction
165            //
166            for(i=0; i<=n-1; i++)
167            {
168                vt = t1[i]/(1+lambda);
169                for(i_=0; i_<=n-1;i_++)
170                {
171                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
172                }
173            }
174        }
175
176
177        /*************************************************************************
178        Inverse matrix update by the Sherman-Morrison formula
179
180        The algorithm updates matrix A^-1 when adding a vector to a column
181        of matrix A.
182
183        Input parameters:
184            InvA        -   inverse of matrix A.
185                            Array whose indexes range within [0..N-1, 0..N-1].
186            N           -   size of matrix A.
187            UpdColumn   -   the column of A whose vector U was added.
188                            0 <= UpdColumn <= N-1
189            U           -   the vector to be added to a column.
190                            Array whose index ranges within [0..N-1].
191
192        Output parameters:
193            InvA        -   inverse of modified matrix A.
194
195          -- ALGLIB --
196             Copyright 2005 by Bochkanov Sergey
197        *************************************************************************/
198        public static void rmatrixinvupdatecolumn(ref double[,] inva,
199            int n,
200            int updcolumn,
201            ref double[] u)
202        {
203            double[] t1 = new double[0];
204            double[] t2 = new double[0];
205            int i = 0;
206            int j = 0;
207            double lambda = 0;
208            double vt = 0;
209            int i_ = 0;
210
211            t1 = new double[n-1+1];
212            t2 = new double[n-1+1];
213           
214            //
215            // T1 = InvA * U
216            // Lambda = v * InvA * U
217            //
218            for(i=0; i<=n-1; i++)
219            {
220                vt = 0.0;
221                for(i_=0; i_<=n-1;i_++)
222                {
223                    vt += inva[i,i_]*u[i_];
224                }
225                t1[i] = vt;
226            }
227            lambda = t1[updcolumn];
228           
229            //
230            // T2 = v*InvA
231            //
232            for(i_=0; i_<=n-1;i_++)
233            {
234                t2[i_] = inva[updcolumn,i_];
235            }
236           
237            //
238            // InvA = InvA - correction
239            //
240            for(i=0; i<=n-1; i++)
241            {
242                vt = t1[i]/(1+lambda);
243                for(i_=0; i_<=n-1;i_++)
244                {
245                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
246                }
247            }
248        }
249
250
251        /*************************************************************************
252        Inverse matrix update by the Sherman-Morrison formula
253
254        The algorithm computes the inverse of matrix A+u*v’ by using the given matrix
255        A^-1 and the vectors u and v.
256
257        Input parameters:
258            InvA    -   inverse of matrix A.
259                        Array whose indexes range within [0..N-1, 0..N-1].
260            N       -   size of matrix A.
261            U       -   the vector modifying the matrix.
262                        Array whose index ranges within [0..N-1].
263            V       -   the vector modifying the matrix.
264                        Array whose index ranges within [0..N-1].
265
266        Output parameters:
267            InvA - inverse of matrix A + u*v'.
268
269          -- ALGLIB --
270             Copyright 2005 by Bochkanov Sergey
271        *************************************************************************/
272        public static void rmatrixinvupdateuv(ref double[,] inva,
273            int n,
274            ref double[] u,
275            ref double[] v)
276        {
277            double[] t1 = new double[0];
278            double[] t2 = new double[0];
279            int i = 0;
280            int j = 0;
281            double lambda = 0;
282            double vt = 0;
283            int i_ = 0;
284
285            t1 = new double[n-1+1];
286            t2 = new double[n-1+1];
287           
288            //
289            // T1 = InvA * U
290            // Lambda = v * T1
291            //
292            for(i=0; i<=n-1; i++)
293            {
294                vt = 0.0;
295                for(i_=0; i_<=n-1;i_++)
296                {
297                    vt += inva[i,i_]*u[i_];
298                }
299                t1[i] = vt;
300            }
301            lambda = 0.0;
302            for(i_=0; i_<=n-1;i_++)
303            {
304                lambda += v[i_]*t1[i_];
305            }
306           
307            //
308            // T2 = v*InvA
309            //
310            for(j=0; j<=n-1; j++)
311            {
312                vt = 0.0;
313                for(i_=0; i_<=n-1;i_++)
314                {
315                    vt += v[i_]*inva[i_,j];
316                }
317                t2[j] = vt;
318            }
319           
320            //
321            // InvA = InvA - correction
322            //
323            for(i=0; i<=n-1; i++)
324            {
325                vt = t1[i]/(1+lambda);
326                for(i_=0; i_<=n-1;i_++)
327                {
328                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
329                }
330            }
331        }
332
333
334        public static void shermanmorrisonsimpleupdate(ref double[,] inva,
335            int n,
336            int updrow,
337            int updcolumn,
338            double updval)
339        {
340            double[] t1 = new double[0];
341            double[] t2 = new double[0];
342            int i = 0;
343            int j = 0;
344            double lambda = 0;
345            double vt = 0;
346            int i_ = 0;
347
348            t1 = new double[n+1];
349            t2 = new double[n+1];
350           
351            //
352            // T1 = InvA * U
353            //
354            for(i_=1; i_<=n;i_++)
355            {
356                t1[i_] = inva[i_,updrow];
357            }
358           
359            //
360            // T2 = v*InvA
361            //
362            for(i_=1; i_<=n;i_++)
363            {
364                t2[i_] = inva[updcolumn,i_];
365            }
366           
367            //
368            // Lambda = v * InvA * U
369            //
370            lambda = updval*inva[updcolumn,updrow];
371           
372            //
373            // InvA = InvA - correction
374            //
375            for(i=1; i<=n; i++)
376            {
377                vt = updval*t1[i];
378                vt = vt/(1+lambda);
379                for(i_=1; i_<=n;i_++)
380                {
381                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
382                }
383            }
384        }
385
386
387        public static void shermanmorrisonupdaterow(ref double[,] inva,
388            int n,
389            int updrow,
390            ref double[] v)
391        {
392            double[] t1 = new double[0];
393            double[] t2 = new double[0];
394            int i = 0;
395            int j = 0;
396            double lambda = 0;
397            double vt = 0;
398            int i_ = 0;
399
400            t1 = new double[n+1];
401            t2 = new double[n+1];
402           
403            //
404            // T1 = InvA * U
405            //
406            for(i_=1; i_<=n;i_++)
407            {
408                t1[i_] = inva[i_,updrow];
409            }
410           
411            //
412            // T2 = v*InvA
413            // Lambda = v * InvA * U
414            //
415            for(j=1; j<=n; j++)
416            {
417                vt = 0.0;
418                for(i_=1; i_<=n;i_++)
419                {
420                    vt += v[i_]*inva[i_,j];
421                }
422                t2[j] = vt;
423            }
424            lambda = t2[updrow];
425           
426            //
427            // InvA = InvA - correction
428            //
429            for(i=1; i<=n; i++)
430            {
431                vt = t1[i]/(1+lambda);
432                for(i_=1; i_<=n;i_++)
433                {
434                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
435                }
436            }
437        }
438
439
440        public static void shermanmorrisonupdatecolumn(ref double[,] inva,
441            int n,
442            int updcolumn,
443            ref double[] u)
444        {
445            double[] t1 = new double[0];
446            double[] t2 = new double[0];
447            int i = 0;
448            int j = 0;
449            double lambda = 0;
450            double vt = 0;
451            int i_ = 0;
452
453            t1 = new double[n+1];
454            t2 = new double[n+1];
455           
456            //
457            // T1 = InvA * U
458            // Lambda = v * InvA * U
459            //
460            for(i=1; i<=n; i++)
461            {
462                vt = 0.0;
463                for(i_=1; i_<=n;i_++)
464                {
465                    vt += inva[i,i_]*u[i_];
466                }
467                t1[i] = vt;
468            }
469            lambda = t1[updcolumn];
470           
471            //
472            // T2 = v*InvA
473            //
474            for(i_=1; i_<=n;i_++)
475            {
476                t2[i_] = inva[updcolumn,i_];
477            }
478           
479            //
480            // InvA = InvA - correction
481            //
482            for(i=1; i<=n; i++)
483            {
484                vt = t1[i]/(1+lambda);
485                for(i_=1; i_<=n;i_++)
486                {
487                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
488                }
489            }
490        }
491
492
493        public static void shermanmorrisonupdateuv(ref double[,] inva,
494            int n,
495            ref double[] u,
496            ref double[] v)
497        {
498            double[] t1 = new double[0];
499            double[] t2 = new double[0];
500            int i = 0;
501            int j = 0;
502            double lambda = 0;
503            double vt = 0;
504            int i_ = 0;
505
506            t1 = new double[n+1];
507            t2 = new double[n+1];
508           
509            //
510            // T1 = InvA * U
511            // Lambda = v * T1
512            //
513            for(i=1; i<=n; i++)
514            {
515                vt = 0.0;
516                for(i_=1; i_<=n;i_++)
517                {
518                    vt += inva[i,i_]*u[i_];
519                }
520                t1[i] = vt;
521            }
522            lambda = 0.0;
523            for(i_=1; i_<=n;i_++)
524            {
525                lambda += v[i_]*t1[i_];
526            }
527           
528            //
529            // T2 = v*InvA
530            //
531            for(j=1; j<=n; j++)
532            {
533                vt = 0.0;
534                for(i_=1; i_<=n;i_++)
535                {
536                    vt += v[i_]*inva[i_,j];
537                }
538                t2[j] = vt;
539            }
540           
541            //
542            // InvA = InvA - correction
543            //
544            for(i=1; i<=n; i++)
545            {
546                vt = t1[i]/(1+lambda);
547                for(i_=1; i_<=n;i_++)
548                {
549                    inva[i,i_] = inva[i,i_] - vt*t2[i_];
550                }
551            }
552        }
553    }
554}
Note: See TracBrowser for help on using the repository browser.