Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Add block-diagonal-matrix construction function #93

Open
mikera opened this Issue · 6 comments

2 participants

@mikera
Owner

Block diagonal matrices are a generally useful type of sparse matrix, and they are supported by various underlying implementations (e.g. Vectorz). See:

It would make sense to have an API function to construct such matrices.

Proposed API:

(block-diagonal-matrix [A B C D ....]

Where A, B, C, D etc. are themselves 2D square matrices which will be arranged along the main diagonal of the resulting matrix.

@prasant94
Collaborator

I have tried to implement this. I added a block-diagonal-matrix function to matrix.clj file and a PBlockDiagonalMatrix protocol to the protocols.clj file.
I am a little confused about the implementation-check method but I just followed the format in the other construction functions and wrote this.

@mikera
Owner

The implementation check is mainly needed to find the current implementation if an implementation is not specified. We need this in matrix construction functions, since we don't always have an existing array to act upon.

@mikera
Owner

I just checked the PR (#104), it looks good in terms of API function and protocol.

We still need:

  • A default implementation (i.e. what code gets called if the implementation doesn't explicitly implement the protocol?) - this needs to construct a block diagonal matrix in a generic way.
  • Tests, ideally as part of the compliance test suite
@prasant94
Collaborator

I have written the default implementation. I am not sure what coerce-param does and if its needed here. The two lines at the end commented out because I didn't know which one to use. Does coerce-param need to be called on the obtained matrix?

(extend-protocol mp/PBlockDiagonalMatrix
  Object
    (block-diagonal-matrix [m blocks]
      (let [num-blocks (mp/dimension-count blocks 0)
          dims-of-block (mp/dimension-count blocks 1)
          dims (* num-blocks dims-of-block)
          zs (vec (repeat dims 0.0))
          dm (vec (for [i (range dims)]
                        (into [] (concat (subvec zs 0 (- i (rem i dims-of-block)))
                                         ((blocks (quot i dims-of-block)) (rem i dims-of-block))
                                         (subvec zs (+ dims-of-block (- i (rem i dims-of-block))))))))]
        ; (mp/coerce-param m dm))))
        ; (dm)))

basically, this

(concat (subvec zs 0 (- i (rem i dims-of-block)))
             ((blocks (quot i dims-of-block)) (rem i dims-of-block))
             (subvec zs (+ dims-of-block (- i (rem i dims-of-block)))))

is concatenating each row of each 2D matrix with the appropriate number of 0s on either side.

@mikera
Owner

This looks close to what is needed.

I think it needs some testing for edge cases: what if the blocks are different sizes, for example? I don't think we can safely assume that dims-of-block is constant.

@prasant94
Collaborator

This one handles those edge cases, but its recursive. Will that be a problem for larger inputs?

(extend-protocol mp/PBlockDiagonalMatrix
  Object
    (block-diagonal-matrix [m blocks]
      (let [aux (fn aux [acc blocks]
                (if (empty? blocks)
                    acc
                    (let [acc-dim (mp/dimension-count acc 0)
                          new-block (blocks 0)
                          new-block-dim (mp/dimension-count new-block 0)
                          new-dim (+ acc-dim new-block-dim)
                          dm (vec (for [i (range new-dim)]
                                       (if (< i acc-dim)
                                           (into [] (concat (acc i)
                                                            (mp/new-vector [] new-block-dim)))
                                           (into [] (concat (mp/new-vector [] acc-dim)
                                                            (new-block (- i acc-dim)))))))]
                          (aux dm (subvec blocks 1)))))]
        (aux [] blocks))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.