matrix/mat64: add `Exp' to TriDense about matrix

[btracey]

The problem is t = NewTriDense(n, computeKind, nil). The pointer itself is passed by value (copied). This doesn't modify what the pointer points to, it changes the pointer altogether. What you really want is *t = *NewTriDense, you want to assign the value the pointer points to for the code to have the intended effect.

[btracey]

I don't understand what's being suggested. The original title is to add Exp to Triangular. That seems fine. By "arithmetic functions" do you mean Add etc? That would be fine too, though it would probably better to have a separate tracking issue.

Is the part about dual numbers an implementation suggestion for Exp, or a suggestion to modify the implementation of TriDense to support the use case of a triangle of dual numbers?

(I do have https://github.com/btracey/hyperdual if you need it)

[ChristopherRabotin]

The suggestion is only to add an Exp function to TriDense. Subsequently, I think it could be interesting to add Sub to TriDense and slightly modify Add of TriDense to simplify the operation if the provided matrix has an even number of rows.

There will be no support for dual numbers, I'll just be using some of their properties as a tool to simplify the number of operations needed.

(P.S.: I think there's a bug in your hyperdual implementation in the Add function. I think line 44 should be a.data[i] += v.)

[ChristopherRabotin]

Sorry for the multiple comments here. I'm encountering an oddity in the implementation, and I suspect that it's due to my total misunderstanding of how reuseAs works. Any hints would be appreciated.

It seems to me that reuseAs prevents modifications of the underlying data, as if there was a deferred function which resets it to the original value.

Here's an output from a basic test that I'm running: the first matrix shown is printed within the new Exp function, and the second one is printed from the main (input matrix is mat64.NewTriDense(2, matrix.Lower, []float64{0, 0, 2, 0})).

⎡1  0⎤
⎣2  1⎦
⎡0  0⎤
⎣0  0⎦

Below is the code executed in this test:

func (t *TriDense) Exp(a Triangular) {
	n, kind := a.Triangle()
	t.reuseAs(a.Triangle())

	if n%2 == 0 {
		// Exponential is trivial. Calculation performed on lower triangular.
		aP := NewTriDense(n, kind, nil)
		isUpper := kind == matrix.Upper
		if isUpper {
			aP.Copy(a.T())
		} else {
			aP.Copy(a)
		}
		computeKind := matrix.Lower
		t = NewTriDense(n, computeKind, nil)
		diag := NewTriDense(n, computeKind, nil)
		eye := NewTriDense(n, computeKind, nil)
		nilDiag := true
		for i := 0; i < n; i++ {
			eye.SetTri(i, i, 1)
			v := aP.At(i, i)
			if v != 0 {
				nilDiag = false
			}
			diag.SetTri(i, i, v)
			aP.SetTri(i, i, 0)
		}
		// TODO(ChristopherRabotin): Use (*TriDense).Add(a, b) once written.
		// Compute exponential which is simply diag*(eye + aP) (if diag != [0]).
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				if i >= j {
					t.SetTri(i, j, eye.At(i, j)+aP.At(i, j))
				}
			}
		}
		if !nilDiag {
			t.MulTri(diag, t)
		}
		if isUpper {
			t.Copy(t.T())
		}
		fmt.Printf("%+v\n", Formatted(t))
	} else {
		// Convert to Dense then convert back to Triangle.
		// The exponential of a triangular matrix is always a triangular matrix
		// of the same kind.
		var triExp Dense
		triExp.Exp(DenseCopyOf(t))
		t.Copy(&triExp) // Store result in receiver.
	}
}

[btracey]

Did you figure out the internal issue (I can describe if not)?

[btracey]

reuseAs doesn't do anything with copying and defer. reuseAs checks the size of the receiver if it is non-zero, and makes it the correct size if it is non-zero.

I'm not exactly sure where you're printing, but t.Copy(t.T()) is wrong. Copy can't be used with itself as a receiver. This is a bug in TriDense.Copy, it should panic like the Dense version does.

[btracey]

As a side note, you should probably working with upper triangular matrices. They're faster in gonum since we use row-major storage.

[ChristopherRabotin]

Okay, I'll switch to using upper triangular and not use copy. I have fixed the internal call issue I originally mentioned (I did a dirty symbolic link but it does the trick). Concerning the main bug, the first matrix printed is the Printf line in Exp, at the end of the scope for even dimensional triangular matrices. The second print is from the actual object on which Exp is called, and after the call. So what happens is that the pointer receiver is never actually modified. I still have yet to understand this.

…

On Fri, Dec 30, 2016, 13:13 Brendan Tracey ***@***.***> wrote: As a side note, you should probably working with upper triangular matrices. They're faster in gonum since we use row-major storage. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#419 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AEma6DivFLel7EVpcPCFwp6GKbSbmiw1ks5rNWXZgaJpZM4LYBf7> .

[kortschak]

So what is the code you have now?

[ChristopherRabotin]

Closed as per #420 (comment).