In Mahout, all of this works pretty well without needing the information
that you are suggesting.
Inversion is not supported, of course, but a solve method on Diagonal can
definitely know what to do. See CholeskyDecomposition for an example.
Addition of a diagonal matrix to a normal matrix works because the abstract
addition method iterates only through non-zeros which gives only the
diagonal elements.
On Fri, Aug 26, 2011 at 4:31 PM, Greg Sterijevski <gsterijev...@gmail.com>wrote:
> Okay, here we go!
>
> Yes, every matrix can be written or considered sparse. My point is that
> knowing the sparsity obeys a pattern (say diagonal) makes your life easier.
>
> Maybe my muddled discussion has added more confusion, an example is in
> order
> (to maximize the entropy, of course). If I was going to build a
> "DiagonalRealMatrix" object,
> I would store the minimum information:
> 1. Data - a vector of n reals
> 2. Indexer -
> public final int index( i, j ){
> if( i == j ) return j;
> return -1;
> }
>
> Presented the information that D is diagonal (through the Java Type System
> or through a type variable that I explicitly set in the constructor), I
> know
> that some operations can be made very simple. Say I want to take the
> inverse
> of D, in place. The matrix function that takes the inverse doesn't need to
> know much. Something like:
>
> new DiagonalInvertFunction() {
> public void eval(double[] data) {
> for( int i = 0 ; i < data.length; i++){
> data[i] = 1.0/data[i];
> }
> }
> }
>
> Would be very quick and scale very well. Contrast this the case where we
> are
> handed a sparse structure and attempt to invert it. Since we are not told
> that it is diagonal, either we need to explicitly test for "diagonality" or
> we proceed with some inversion method that can handle arbitrary forms of
> sparseness. Both approaches (testing or using a sparse invertor) must be
> more expensive? No?
>
> Say that we are adding the diagonal above to a generic matrix (this is a
> bit
> outside of the scope of the discussion). I know, because of the pattern of
> the matrix, I need only make n additions, instead of looping naively
> through
> based on some striding scheme.
>
> Where my idea falls apart is precisely in operations like the one you cite
> below. Say that the user wants to calculate, elementwise, the cosine of
> each
> element in this diagonal matrix. Here the situation is messy because the
> pattern of the data changes. The tidy vector which stores just the diagonal
> elements, is not appropriate since the diagonal matrix now becomes a
> general
> matrix.
>
> In the end, I think I see your point, but, its an empirical question as to
> whether my examples (as archetypes) are more likely to occur in practice
> versus your example.
>
> -Greg
>
>
> On Fri, Aug 26, 2011 at 1:41 PM, Ted Dunning <ted.dunn...@gmail.com>
> wrote:
>
> > On Fri, Aug 26, 2011 at 7:38 AM, Greg Sterijevski <
> gsterijev...@gmail.com
> > >wrote:
> >
> > > Ted,
> > >
> > > When you say
> > >
> > > "Functions are good, but giving a tiny bit
> > > more information to the function is also a great idea"
> > >
> > > do you mean information on indexing and shape of the data?
> > >
> >
> > I meant the location of the element. For instance, it would be super
> > simple
> > to define Hilbert's matrix if you had extra arguments to the function:
> >
> > Matrix a = new DenseMatrix(10,10).assign(new ElementFunction() {
> > public double eval(int i, int j, double value) { return 1.0 / (1 +
> i
> > + j); }
> > }
> >
> > Without the indexes, you are reduced to writing loops. With the indexes,
> > you have a very concise functional definition.
> >
> > I have missed this more than once in Mahout.
> >
> >
> >
> > > One thought I had, I am not sure if this is 100% applicable is the
> > > following:
> > >
> > > 1. You have two types of data (in matrix) : structured and sparse.
> > > 2. Structured data is a collection of elements that fits a certain
> > pattern.
> > > Some of which are
> > > i.) General matrix, stored rowwise
> > > ii.) General matrix, stored columnwise
> > > iii.) Diagonal matrix
> > >
> >
> > Note that this is really just a special kind of sparse matrix.
> >
> >
> > > iii.) Upper/Lower triangular
> > > iv.) Symmetric stored in either upper triangular or lower triangular
> > > compressed format
> > > v.) Banded
> > >
> >
> > And so is this a kind of sparse matrix.
> >
> >
> > > The visitor should know the pattern and be given an indexer function.
> >
> >
> > Why?
> >
> > Can you adduce a use case for knowing the pattern inside the visitor?
> > Isn't
> > it true that most special cases fall into two situations:
> >
> > - the caller knows and can instantiate a special function fit for the
> > purpose? For instance, in computing an SVD, you *know* that you have a
> > bi-diagonal form at some point. You don't need to discover this.
> >
> > - standard sparse techniques work fine. For instance, if you do have a
> > banded matrix that just does the normal sparse thing for adding and
> > multiplying, don't you get what you want?
> >
> > While I think that this breakdown pretty much covers it for element-wise
> > functions, this clearly doesn't apply to functions of an entire matrix.
> If
> > you are doing an SVD, you may find it very advantageous to handle, say
> > orthornormal matrices or diagonal offsets of a previous decomposition or
> > banded forms specially.
> >
> >
> >
> > > The pattern designation eliminates superfluous calls to the indexer.
> You
> > > would
> > > not ask for element i,j when i ne j, if the matrix is diagonal.
> >
> >
> > I think that I need to see a use case to respond.
> >
> > In general, you *do* need to see all the elements. For instance, suppose
> > that I am setting all elements of a matrix to (a_ij + 1)^2 or taking the
> > cosine of all elements.
> >
> > It might be that you want to have a specialized assignment for just some
> > elements. Or you might want to have a strange kind of view that just
> > doesn't have some elements.
> >
> >
> > > ... The sparse structure is very important, but since you are unsure of
> > how
> > > the
> > > sparseness occurs, whether there is a pattern or not, you could make
> > blind
> > > calls to the indexer, or the indexer itself could be class (in this
> case)
> > > which returns to you a collection of acceptable tuples ( eg, element
> > > locations which might or might not be zero...)
> > >
> >
> > This is abstract enough that I don't know how to relate this to use
> cases.
> > Can you provide a specific example?
> >
> > If I want to do something to just the diagonal elements, I generally use
> > the
> > viewDiagonal method. I can imagine a viewNonZero(double epsilon) which
> > gives a view that is a matrix whose purpose in life is to restrict the
> > scope
> > of the assign operation, but I usually just use the iterateNonZero
> element
> > iterator for those cases.
> >
> >
> > > In cases of very small problems, Luc's 3x3 or 6x6 matrix, it makes
> sense
> > to
> > > have subclasses of the general matrix, with rows and columns fixed to
> be
> > 3
> > > or 6, respectively.
> > >
> >
> > Absolutely true. It might even pay to automatically generate special
> cases
> > for known important sizes (3x3, yes, 4x4, yes, 6x6, yes, 3x5, probably
> > not).
> >
>
