SparseMatrix


1   /*
2    * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
3    * Copyright (C) 2006 - JScience (http://jscience.org/)
4    * All rights reserved.
5    * 
6    * Permission to use, copy, modify, and distribute this software is
7    * freely granted, provided that this notice is preserved.
8    */
9   package org.jscience.mathematics.vectors;
10  
11  import java.util.Iterator  ;
12  import java.util.List  ;
13  
14  import javolution.lang.MathLib;
15  import javolution.util.FastComparator;
16  import javolution.util.FastMap;
17  import javolution.util.FastTable;
18  import javolution.util.Index;
19  
20  import org.jscience.mathematics.structures.Field;
21  
22  /**
23   * <p> This class represents a matrix made of {@link SparseVector sparse
24   *     vectors} (as rows). To create a sparse matrix made of column vectors the 
25   *     {@link #transpose} method can be used. 
26   *     For example:[code]
27   *        SparseVector<Rational> column0 = SparseVector.valueOf(...);
28   *        SparseVector<Rational> column1 = SparseVector.valueOf(...);
29   *        SparseMatrix<Rational> M = SparseMatrix.valueOf(column0, column1).transpose();
30   *     [/code]</p>
31   * <p> As for any concrete {@link org.jscience.mathematics.structures.Structure
32   *     structure}, this class is declared <code>final</code> (otherwise most
33   *     operations would have to be overridden to return the appropriate type).
34   *     Specialized dense matrix should sub-class {@link Matrix} directly.
35   *     For example:[code]
36   *        // Extension through composition.
37   *        final class BandMatrix <F extends Field<F>> extends Matrix<F> {
38   *             private SparseMatrix<F> _value;
39   *             ...
40   *             public BandMatrix opposite() { // Returns the right type.
41   *                 return BandMatrix.valueOf(_value.opposite());
42   *             }
43   *             ...
44   *        }[/code]
45   *     </p>   
46   * @author <a HREF="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
47   * @version 3.3, January 2, 2007
48   */
49  public final class SparseMatrix<F extends Field<F>> extends Matrix<F> {
50  
51      /**
52       * Holds the number of columns n or the number of rows m if transposed.
53       */
54      int _n;
55  
56      /**
57       * Holds the zero.
58       */
59      F _zero;;
60  
61      /**
62       * Indicates if this matrix is transposed (the rows are then the columns).
63       */
64      boolean _transposed;
65  
66      /**
67       * Holds this matrix rows (or columns when transposed).
68       */
69      final FastTable<SparseVector<F>> _rows = new FastTable<SparseVector<F>>();
70  
71      /**
72       * Returns the sparse square matrix having the specified diagonal
73       * vector. This method is typically used to create an identity matrix.
74       * For example:[code]
75       *      SparseMatrix<Real> IDENTITY = Matrix.valueOf(
76       *           DenseVector.valueOf({Real.ONE, Real.ONE, Real.ONE}), Real.ZERO);
77       * [/code]          
78       *
79       * @param  diagonal the diagonal vector.
80       * @param  zero value of non-diagonal elements.
81       * @return a square matrix with <code>diagonal</code> on the diagonal and
82       *         <code>zero</code> elsewhere.
83       */
84      public static <F extends Field<F>> SparseMatrix<F> valueOf(
85              Vector<F> diagonal, F zero) {
86          int n = diagonal.getDimension();
87          SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
88          for (int i = 0; i < n; i++) {
89              SparseVector<F> row = SparseVector.valueOf(n, zero, i, diagonal
90                      .get(i));
91              M._rows.add(row);
92          }
93          return M;
94      }
95  
96      /**
97       * Returns a sparse matrix holding the specified row vectors 
98       * (column vectors if {@link #transpose transposed}).
99       *
100      * @param rows the row vectors.
101      * @return the matrix having the specified rows.
102      * @throws DimensionException if the rows do not have the same dimension.
103      */
104     public static <F extends Field<F>> SparseMatrix<F> valueOf(
105             SparseVector<F>... rows) {
106         final int n = rows[0]._dimension;
107         final F zero = rows[0]._zero;
108         SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
109         for (int i = 0, m = rows.length; i < m; i++) {
110             SparseVector<F> rowi = rows[i];
111             if (rowi._dimension != n)
112                 throw new DimensionException(
113                         "All vectors must have the same dimension.");
114             if (!zero.equals(rowi._zero))
115                 throw new DimensionException(
116                         "All vectors must have the same zero element.");
117             M._rows.add(rowi);
118         }
119         return M;
120     }
121 
122     /**
123      * Returns a sparse matrix holding the row vectors from the specified 
124      * collection (column vectors if {@link #transpose transposed}).
125      *
126      * @param rows the list of row vectors.
127      * @return the matrix having the specified rows.
128      * @throws DimensionException if the rows do not have the same dimension.
129      */
130     public static <F extends Field<F>> SparseMatrix<F> valueOf(
131             List  <SparseVector<F>> rows) {
132         final int n = rows.get(0)._dimension;
133         final F zero = rows.get(0)._zero;
134         SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
135         Iterator  <SparseVector<F>> iterator = rows.iterator();
136         for (int i = 0, m = rows.size(); i < m; i++) {
137             SparseVector<F> rowi = iterator.next();
138             if (rowi.getDimension() != n)
139                 throw new DimensionException(
140                         "All vectors must have the same dimension.");
141             if (!zero.equals(rowi._zero))
142                 throw new DimensionException(
143                         "All vectors must have the same zero element.");
144             M._rows.add(rowi);
145         }
146         return M;
147     }
148 
149     /**
150      * Returns a sparse matrix equivalent to the specified matrix but with 
151      * the zero elements removed using the default object equality comparator.
152      *
153      * @param that the matrix to convert.
154      * @param zero the zero element for the sparse vector to return.
155      * @return <code>SparseMatrix.valueOf(that, zero, FastComparator.DEFAULT)</code> or a dense matrix holding the same elements
156      */
157     public static <F extends Field<F>> SparseMatrix<F> valueOf(Matrix<F> that, F zero) {
158         return SparseMatrix.valueOf(that, zero, FastComparator.DEFAULT);
159     }
160 
161     /**
162      * Returns a sparse matrix equivalent to the specified matrix but with 
163      * the zero elements removed using the specified object equality comparator.
164      *
165      * @param that the matrix to convert.
166      * @param zero the zero element for the sparse vector to return.
167      * @param comparator the comparator used to determinate zero equality. 
168      * @return <code>that</code> or a dense matrix holding the same elements
169      *         as the specified matrix.
170      */
171     public static <F extends Field<F>> SparseMatrix<F> valueOf(Matrix<F> that,
172             F zero, FastComparator<? super F> comparator) {
173         if (that instanceof SparseMatrix)
174             return (SparseMatrix<F>) that;
175         int n = that.getNumberOfColumns();
176         int m = that.getNumberOfRows();
177         SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
178         for (int i = 0; i < m; i++) {
179             SparseVector<F> rowi = SparseVector.valueOf(that.getRow(i), zero,
180                     comparator);
181             M._rows.add(rowi);
182         }
183         return M;
184     }
185 
186     /**
187      * Returns the value of the non-set elements for this sparse matrix.
188      * 
189      * @return the element corresponding to zero.
190      */
191     public F getZero() {
192         return _zero;
193     }
194 
195     @Override  
196     public int getNumberOfRows() {
197         return _transposed ? _n : _rows.size();
198     }
199 
200     @Override  
201     public int getNumberOfColumns() {
202         return _transposed ? _rows.size() : _n;
203     }
204 
205     @Override  
206     public F get(int i, int j) {
207         return _transposed ? _rows.get(j).get(i) : _rows.get(i).get(j);
208     }
209 
210     @Override  
211     public SparseVector<F> getRow(int i) {
212         if (!_transposed)
213             return _rows.get(i);
214         // Else transposed.
215         int n = _rows.size();
216         int m = _n;
217         if ((i < 0) || (i >= m))
218             throw new DimensionException();
219         SparseVector<F> V = SparseVector.newInstance(n, _zero);
220         for (int j = 0; j < n; j++) {
221             SparseVector<F> row = _rows.get(j);
222             F e = row._elements.get(Index.valueOf(i));
223             if (e != null) {
224                 V._elements.put(Index.valueOf(j), e);
225             }
226         }
227         return V;
228     }
229 
230     @Override  
231     public SparseVector<F> getColumn(int j) {
232         if (_transposed)
233             return _rows.get(j);
234         int m = _rows.size();
235         if ((j < 0) || (j >= _n))
236             throw new DimensionException();
237         SparseVector<F> V = SparseVector.newInstance(_n, _zero);
238         for (int i = 0; i < m; i++) {
239             SparseVector<F> row = _rows.get(i);
240             F e = row._elements.get(Index.valueOf(j));
241             if (e != null) {
242                 V._elements.put(Index.valueOf(i), e);
243             }
244         }
245         return V;
246     }
247 
248     @Override  
249     public SparseVector<F> getDiagonal() {
250         int m = this.getNumberOfRows();
251         int n = this.getNumberOfColumns();
252         int dimension = MathLib.min(m, n);
253         SparseVector<F> V = SparseVector.newInstance(_n, _zero);
254         for (int i = 0; i < dimension; i++) {
255             SparseVector<F> row = _rows.get(i);
256             F e = row._elements.get(Index.valueOf(i));
257             if (e != null) {
258                 V._elements.put(Index.valueOf(i), e);
259             }
260         }
261         return V;
262     }
263 
264     @Override  
265     public SparseMatrix<F> opposite() {
266         SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed);
267         for (int i = 0, p = _rows.size(); i < p; i++) {
268             M._rows.add(_rows.get(i).opposite());
269         }
270         return M;
271     }
272 
273     @Override  
274     public SparseMatrix<F> plus(Matrix<F> that) {
275         if (this.getNumberOfRows() != that.getNumberOfRows())
276             throw new DimensionException();
277         SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed);
278         for (int i = 0, p = _rows.size(); i < p; i++) {
279             M._rows.add(_rows.get(i).plus(
280                     _transposed ? that.getColumn(i) : that.getRow(i)));
281         }
282         return M;
283     }
284 
285     @Override  
286     public SparseMatrix<F> minus(Matrix<F> that) { // Returns more specialized type.
287         return this.plus(that.opposite());
288     }
289 
290     @Override  
291     public SparseMatrix<F> times(F k) {
292         SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed);
293         for (int i = 0, p = _rows.size(); i < p; i++) {
294             M._rows.add(_rows.get(i).times(k));
295         }
296         return M;
297     }
298 
299     @Override  
300     public SparseVector<F> times(Vector<F> v) {
301         if (v.getDimension() != this.getNumberOfColumns())
302             throw new DimensionException();
303         final int m = this.getNumberOfRows();
304         SparseVector<F> V = SparseVector.newInstance(m, _zero);
305         for (int i = 0; i < m; i++) {
306             F e = this.getRow(i).times(v);
307             if (!_zero.equals(e)) {
308                 V._elements.put(Index.valueOf(i), e);
309             }
310         }
311         return V;
312     }
313 
314     @Override  
315     public SparseMatrix<F> times(Matrix<F> that) {
316         final int m = this.getNumberOfRows();
317         final int n = this.getNumberOfColumns();
318         final int p = that.getNumberOfColumns();
319         if (that.getNumberOfRows() != n)
320             throw new DimensionException();
321         // Creates a mxp matrix in transposed form (p columns vectors of size m)
322         SparseMatrix<F> M = SparseMatrix.newInstance(m, _zero, true);   
323         for (int j = 0; j < p; j++) { 
324             Vector<F> thatColj = that.getColumn(j);
325             SparseVector<F> column = SparseVector.newInstance(m, _zero);
326             M._rows.add(column); // M is transposed.
327             for (int i = 0; i < m; i++) {
328                 SparseVector<F> thisRowi = this.getRow(i);
329                 F e = thisRowi.times(thatColj);
330                 if (!_zero.equals(e)) {
331                     column._elements.put(Index.valueOf(i), e);
332                 }
333             }
334         }
335         return M;
336     }
337 
338     @Override  
339     public SparseMatrix<F> inverse() {
340         if (!isSquare())
341             throw new DimensionException("Matrix not square");
342         F detInv = this.determinant().inverse();
343         SparseMatrix<F> A = this.adjoint();
344         // Multiply adjoint elements with 1 / determinant.
345         for (int i = 0, m = A._rows.size(); i < m; i++) {
346             SparseVector<F> row = A._rows.get(i);
347             for (FastMap.Entry<Index, F> e = row._elements.head(), end = row._elements
348                     .tail(); (e = e.getNext()) != end;) {
349                 F element = e.getValue();
350                 e.setValue(detInv.times(element));
351             }
352         }
353         return A;
354     }
355 
356     @Override  
357     public F determinant() {
358         if (!isSquare())
359             throw new DimensionException("Matrix not square");
360         if (_n == 1)
361             return this.get(0, 0);
362         // Expansion by minors (also known as Laplacian)
363         // This algorithm is division free but too slow for dense matrix.
364         SparseVector<F> row0 = this.getRow(0);
365         F det = null;
366         for (FastMap.Entry<Index, F> e = row0._elements.head(), end = row0._elements
367                 .tail(); (e = e.getNext()) != end;) {
368             int i = e.getKey().intValue();
369             F d = e.getValue().times(cofactor(0, i));
370             if (i % 2 != 0) {
371                 d = d.opposite();
372             }
373             det = (det == null) ? d : det.plus(d);
374         }
375         return det == null ? _zero : det;
376     }
377 
378     @Override  
379     public Matrix<F> solve(Matrix<F> y) {
380         return this.inverse().times(y);
381     }
382 
383     @Override  
384     public SparseMatrix<F> transpose() {
385         SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, !_transposed);
386         M._rows.addAll(this._rows);
387         return M;
388     }
389 
390     @Override  
391     public F cofactor(int i, int j) {
392         if (_transposed) {
393             int k = i;
394             i = j;
395             j = k; // Swaps i,j
396         }
397         int m = _rows.size();
398         SparseMatrix<F> M = SparseMatrix.newInstance(m - 1, _zero, _transposed);
399         for (int k1 = 0; k1 < m; k1++) {
400             if (k1 == i)
401                 continue;
402             SparseVector<F> row = _rows.get(k1);
403             SparseVector<F> V = SparseVector.newInstance(_n - 1, _zero);
404             M._rows.add(V);
405             for (FastMap.Entry<Index, F> e = row._elements.head(), end = row._elements
406                     .tail(); (e = e.getNext()) != end;) {
407                 int index = e.getKey().intValue();
408                 if (index < j) {
409                     V._elements.put(e.getKey(), e.getValue());
410                 } else if (index > j) { // Position shifted (index minus one).
411                     V._elements.put(Index.valueOf(index - 1), e.getValue());
412                 } // Else don't copy element at j.
413             }
414         }
415         return M.determinant();
416     }
417 
418     @Override  
419     public SparseMatrix<F> adjoint() {
420         SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed);
421         int m = _rows.size();
422         for (int i = 0; i < m; i++) {
423             SparseVector<F> row = SparseVector.newInstance(_n, _zero);
424             M._rows.add(row);
425             for (int j = 0; j < _n; j++) {
426                 F cofactor = _transposed ? cofactor(j, i) : cofactor(i, j);
427                 if (!_zero.equals(cofactor)) {
428                     row._elements
429                             .put(Index.valueOf(j),
430                                     ((i + j) % 2 == 0) ? cofactor : cofactor
431                                             .opposite());
432                 }
433             }
434         }
435         return M.transpose();
436     }
437 
438     @Override  
439     public SparseMatrix<F> tensor(Matrix<F> that) {
440         final int thism = this.getNumberOfRows();
441         final int thisn = this.getNumberOfColumns();
442         final int thatm = that.getNumberOfRows();
443         final int thatn = that.getNumberOfColumns();
444         int n = thisn * thatn; // Number of columns,
445         int m = thism * thatm; // Number of rows.
446         SparseMatrix<F> M = SparseMatrix.newInstance(n, _zero, false);
447         for (int i=0; i < m; i++) { // Row index.
448             final int i_rem_thatm = i % thatm;
449             final int i_div_thatm = i / thatm;
450             SparseVector<F> row = SparseVector.newInstance(n, _zero);
451             M._rows.add(row);
452             SparseVector<F> thisRow = this.getRow(i_div_thatm);
453             for (FastMap.Entry<Index, F> e = thisRow._elements.head(),
454                     end = thisRow._elements.tail(); (e = e.getNext()) != end;) {
455                 F a = e.getValue();
456                 int j = e.getKey().intValue();
457                 for (int k=0; k < thatn; k++) {
458                     F b = that.get(i_rem_thatm, k);
459                     if (!b.equals(_zero)) {
460                         row._elements.put(Index.valueOf(j * thatn + k), a.times(b));
461                     }
462                 }
463             }
464         }
465         return M;
466     }
467 
468     @Override  
469     public SparseVector<F> vectorization() {
470         SparseVector<F> V = SparseVector.newInstance(_n
471                 * this.getNumberOfRows(), _zero);
472         int offset = 0;
473         for (int j = 0, n = this.getNumberOfColumns(); j < n; j++) {
474             SparseVector<F> column = this.getColumn(j);
475             for (FastMap.Entry<Index, F> e = column._elements.head(), end = column._elements
476                     .tail(); (e = e.getNext()) != end;) {
477                 V._elements.put(Index.valueOf(e.getKey().intValue() + offset),
478                         e.getValue());
479             }
480             offset += this.getNumberOfRows();
481         }
482         return V;
483     }
484 
485     @Override  
486     public boolean move(ObjectSpace os) {
487         if (super.move(os)) {
488             // The rows table itself is intrinsic (all final members are)
489             // and implicitely moves with this matrix.
490             // But it may refer to locally/stack allocated vectors which
491             // need to be moved along.
492             for (int i = 0, n = _rows.size(); i < n; i++) {
493                 _rows.get(i).move(os);
494             }
495             _zero.move(os);
496             return true;
497         }
498         return false;
499     }
500 
501     ///////////////////////
502     // Factory creation. //
503     ///////////////////////
504 
505     @SuppressWarnings  ("unchecked")
506     static <F extends Field<F>> SparseMatrix<F> newInstance(int n, F zero,
507             boolean transposed) {
508         SparseMatrix<F> M = FACTORY.object();
509         M._n = n;
510         M._zero = zero;
511         M._transposed = transposed;
512         return M;
513     }
514 
515     private static Factory<SparseMatrix> FACTORY = new Factory<SparseMatrix>() {
516         @Override  
517         protected SparseMatrix create() {
518             return new SparseMatrix();
519         }
520 
521         @Override  
522         protected void cleanup(SparseMatrix matrix) {
523             matrix._rows.reset();
524         }
525     };
526 
527     private SparseMatrix() {
528     }
529 
530     private static final long serialVersionUID = 1L;
531 
532 }
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags