mat: add basis for band matrix type #95

kortschak · 2017-06-21T03:57:30Z

This adds a starting point for building. I'm not currently happy with matrix construction since it's possible as the code stands to change the test case from 6x6 to 6x4 and no other changes without failure. I don't see a way to address this.

@btracey Please take a look.

btracey · 2017-06-23T16:53:17Z

mat/band.go

+// data is used as the backing slice, and changes to the elements of the returned BandDense
+// will be reflected in data. If neither of these is true, NewBandDense will panic.
+//
+// The data must be arranged in row-major order. Only the values in the band portion


I think this needs a definition of "row-major order", as if you haven't seen banded matrices it's not obvious what the layout is. Something like a copy and paste from the BLAS documentation.

btracey · 2017-06-23T16:55:56Z

mat/band.go

+	mat blas64.Band
+}
+
+type Band interface {


I think we want the interface to be Banded rather than Band to match Triangular and Symmetric. I'm not sure if we want the type to be BandedDense instead of BandDense

btracey · 2017-06-23T16:57:48Z

mat/band_test.go

+}
+
+func TestBandAtSet(t *testing.T) {
+	band := NewBandDense(6, 6, 1, 2, []float64{


It would be nice to have a test that just accesses all of the matrix elements to make sure the index tests are correct.

btracey · 2017-06-23T17:06:01Z

I also do not see a way to address the allocation problem.

Two general concerns moving forward (not a problem with this PR).

Do we have / could we have a way of maintaining a list of the special cases that we have implemented? It would be nice to have some big table where we could see gaps that need to be filled.
Here is a nice property of banded matrices. The Cholesky decomposition of a (symmetric) banded matrix is also banded. It would be great if the Cholesky type could deal with that and "just work". In particular, somehow making, say, distmv.Normal agnostic to the representation. I believe the constraint on that particular issue is having the "Lower" extraction be agnostic.

kortschak · 2017-06-30T01:58:24Z

PTAL

vladimir-ch · 2017-07-01T06:57:41Z

mat/band.go

+// TransposeBand is a type for performing an implicit transpose of a band
+// matrix. It implements the Banded interface, returning values from the
+// transpose of the matrix within.
+type TransposeBand struct {


Should this be TransposeBanded? In general, from a quick look, I don't understand the rules when Band or Banded is used.

I did write it like that, but reverted because of length. Various people call these matrices band or banded. Here Banded is the interface and Band is the theoretical notion. I don't really mind either way.

btracey · 2017-07-06T04:19:55Z

mat/band.go

+	Matrix
+	// Triangular returns the number of rows/columns in the matrix and its
+	// orientation.
+	Bandwidth() (kl, ku int)


Should this be plural?

The band's width is kl+ku+1, so no, I don't think so.

But the doc comment should not say Triangular!

btracey · 2017-07-06T04:21:36Z

mat/band.go

+
+	// TBanded is the equivalent of the T() method in the Matrix interface but
+	// guarantees the transpose is of banded type.
+	TBanded() Banded


This should be TBand to match TTri and the rest of the "abbreviation" conventions.

btracey · 2017-07-06T04:27:04Z

mat/band_test.go

+		if cols != test.c {
+			t.Errorf("unexpected number of cols for test %d: got: %d want: %d", i, cols, test.c)
+		}
+		if !reflect.DeepEqual(band, test.mat) {


Add a test here with the *Dense matrix and also use mat.Equal

Which *Dense matrix?

I would write out the banded matrix as a *Dense with the zeros, and then test that equality.

Done, I think.

I was hoping for

matDense: NewDense(6, 6, []float64{ 1, 2, 3, 0, 0, 0, 4, 5, 6, 7, 0, 0, 0, 8, 9, 10, 11, 0, 0, 0, 12, 13, 14, 15, 0, 0, 0, 16, 17, 18, 0, 0, 0, 0, 19, 20, })

and then comparing that with test.band

Ah, OK. Done.

btracey · 2017-07-06T04:28:27Z

mat/band_test.go

+	}
+}
+
+// bandImplicit is an implicit band matrix returning val(i, j)


Could you add these types to list_test as well?

Can I leave that for you? That system scares me a little.

Me too. Sure.

kortschak · 2017-07-06T04:49:10Z

I will lose approval with these changes, so I'll rebase on to master while I'm here.

kortschak · 2017-07-06T05:00:24Z

PTAL

btracey · 2017-07-06T05:02:46Z

It should also probably actually say what the function does. Sorry for missing that.

…

On Jul 5, 2017, at 11:01 PM, Dan Kortschak ***@***.***> wrote: @kortschak commented on this pull request. In mat/band.go <#95 (comment)>: > + bandDense *BandDense + _ Matrix = bandDense + _ Banded = bandDense + _ RawBand = bandDense +) + +// BandDense represents a band matrix in dense storage format. +type BandDense struct { + mat blas64.Band +} + +type Banded interface { + Matrix + // Triangular returns the number of rows/columns in the matrix and its + // orientation. + Bandwidth() (kl, ku int) But the doc comment should not say Triangular! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#95 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADgqWxl9RYjs2UJULoi_wib5qesIcXbSks5sLGoRgaJpZM4OAbXx>.

kortschak · 2017-07-06T05:04:59Z

OK, that and another missed typo are fixed. PTAL

kortschak · 2017-07-06T05:42:51Z

I'll make those changes later today. This will invalidate the approval.

kortschak · 2017-07-06T07:44:01Z

PTAL

vladimir-ch · 2017-07-06T15:11:31Z

mat/band.go

+//    12 13 14 15
+//    16 17 18  *
+//    19 20  *  *
+// which is given to the passed to NewBandDense as []float64{*, 1, 2, 3, 4, ...}.


which is passed to NewBandDense ... ?
Also, wouldn't it be worth to write that in this example kl == 1 and ku == 2?

vladimir-ch · 2017-07-06T15:12:54Z

mat/band.go

+	return TransposeBand{b}
+}
+
+// RawMatrix returns the underlying blas64.Band used by the receiver.


vladimir-ch · 2017-07-06T15:13:01Z

mat/band.go

+	return b.mat.Rows, b.mat.Cols
+}
+
+// Dims returns the upper and lower bandwidths of the matrix.


// Bandwidth returns the lower and upper bandwidths ...

vladimir-ch · 2017-07-06T15:18:27Z

mat/band_test.go

+			}
+		}
+	}
+	// Do that same thing via a vall to Equal.


kortschak · 2017-07-06T22:02:51Z

PTAL

kortschak · 2017-07-07T23:35:53Z

ping @btracey since you are up.

kortschak · 2017-07-08T01:51:53Z

mat/band.go

+// a new slice is allocated for the backing slice. If len(data) == min(r, c+kl)*(kl+ku+1),
+// data is used as the backing slice, and changes to the elements of the returned
+// BandDense will be reflected in data. If neither of these is true, NewBandDense
+// will panic. kl must be at least zero and less r, and ku must be at least zero and


Note comments on kl and ku now.

Does this behavior match Lapack? I think we want to be compliant here. I've only started looking at lapack banded implementations, but it would be nice to not have to play too many games (that said, we are stricter for Vector and it's worked out fine)

The documentation says that the * elements are not referenced. If this causes a failure, I think then that the LAPACK implementation is buggy.

Looking at how we do this, I suspect it will fail for the gonum implementation, for example here. However, I think that is a result of our slice optimisations; the Fortran should not fail.

Thanks. I agree. We/I should probably revisit the banded BLAS implementations before we're 6 layers deep in Lapack trying to figure out what's going on.

There are two ways to deal with those, either continue if the proposed length is 0, or probably preferably, make the same assessment that I am doing here.

btracey · 2017-07-08T05:18:57Z

mat/band.go

@@ -136,6 +136,12 @@ func NewBandDense(r, c, kl, ku int, data []float64) *BandDense {
 	}
 }

+// NewDiagonal is a convenience function that returns a diagonal matrix represented by a
+// BandDense. The length of data must be min(r, c) otherwise NewDiagonal will panic.
+func NewDiagonal(r, c int, data []float64) *BandDense {


Please leave this out of this pull request.

I think we want a SymBanded interface / type. It makes a huge difference in distmv.NewNormal (O(n^3) to O(n) for computing the Cholesky), and is probably a common case in the spatial statistics world also. Typically, diagonal matrices are also square, so I think we would want this to return a SymBandDense by default.

The sigma in SVD is potentially non-square. Similarly in GSVD. NewDiagonalBand and NewDiagonalSym? The reason this is here is purely for legibility to draw attention to the fact that the result is a diagonal.

Yes, I agree with the intent. Wikipedia suggests that the term Diagonal "usually refers to square matrices", and calls non-square ones "rectangular diagonal matrices".

I would go with NewDiagonal for the symmetric case, and NewDiagonalBand or NewDiagonalRect for the rectangular case.

data should be allowed to be nil, just like the other NewXxx functions. I commonly create a new matrix then use a loop to set the diagonal elements (and fortunately these will be much more efficient soon).

Sorry, I'm an idiot. It can already be nil.

An alternative is to have NewDiagonal(r, c int, data []float64) Matrix and return either a *BandDense or a *SymBandDense depending on whether r == c.

I'd rather have the two functions that return their specific types and not force type assertions.

We could also have NewDiagSq and NewDiagRect

I'll rename this to NewDiagonalRect.

btracey reviewed Jun 23, 2017

View reviewed changes

kortschak force-pushed the banded branch from c014515 to c4e05c8 Compare June 30, 2017 01:57

vladimir-ch reviewed Jul 1, 2017

View reviewed changes

btracey previously approved these changes Jul 6, 2017

View reviewed changes

kortschak dismissed btracey’s stale review via f3a39dc July 6, 2017 05:00

kortschak force-pushed the banded branch from c4e05c8 to f3a39dc Compare July 6, 2017 05:00

kortschak force-pushed the banded branch 2 times, most recently from 0c50f60 to 5872b1e Compare July 6, 2017 05:04

btracey previously approved these changes Jul 6, 2017

View reviewed changes

kortschak dismissed btracey’s stale review via 8e0e5a3 July 6, 2017 07:43

kortschak force-pushed the banded branch from 5872b1e to 8e0e5a3 Compare July 6, 2017 07:43

btracey previously approved these changes Jul 6, 2017

View reviewed changes

vladimir-ch reviewed Jul 6, 2017

View reviewed changes

kortschak dismissed btracey’s stale review via ffc4d79 July 6, 2017 22:02

kortschak force-pushed the banded branch from 8e0e5a3 to ffc4d79 Compare July 6, 2017 22:02

kortschak mentioned this pull request Jul 7, 2017

stat/spatial: use mat.RowNonZeroDoer #124

Merged

kortschak force-pushed the banded branch from 18da911 to 86c7bdd Compare July 8, 2017 01:44

kortschak added 2 commits July 8, 2017 11:20

mat: add basis for band matrix type

b92b873

mat: tighten bound on band data

0ef4235

kortschak force-pushed the banded branch from 86c7bdd to 87d408a Compare July 8, 2017 01:51

kortschak commented Jul 8, 2017

View reviewed changes

kortschak mentioned this pull request Jul 8, 2017

mat: add Doer interface and implementations #126

Merged

btracey reviewed Jul 8, 2017

View reviewed changes

mat: add NewDiagonalRect convenience function

a6d9d8f

kortschak force-pushed the banded branch from 87d408a to a6d9d8f Compare July 8, 2017 08:32

btracey approved these changes Jul 8, 2017

View reviewed changes

kortschak merged commit c08ba23 into master Jul 8, 2017

kortschak deleted the banded branch July 8, 2017 22:13

mat: add basis for band matrix type #95

mat: add basis for band matrix type #95

Conversation

kortschak commented Jun 21, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

btracey commented Jun 23, 2017

kortschak commented Jun 30, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak commented Jul 6, 2017

kortschak commented Jul 6, 2017

btracey commented Jul 6, 2017 via email

kortschak commented Jul 6, 2017

kortschak commented Jul 6, 2017 via email

kortschak commented Jul 6, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak commented Jul 6, 2017

kortschak commented Jul 7, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment