Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/inverseupdate.cs @ 3932

Last change on this file since 3932 was 3839, checked in by mkommend, 15 years ago

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