scripts/nn/layers/conv2d.dml

#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

/*
 * 2D Convolutional layer.
 */
source("scripts/nn/util.dml") as util

forward = function(matrix[double] X, matrix[double] W, matrix[double] b,
                   int C, int Hin, int Win, int Hf, int Wf,
                   int strideh, int stridew, int padh, int padw)
    return (matrix[double] out, int Hout, int Wout) {
  /*
   * Computes the forward pass for a 2D spatial convolutional layer with
   * F filters.  The input data has N examples, each represented as a 3D
   * volume unrolled into a single vector.
   *
   * This implementation uses `im2col` internally for each image to
   * extract local image regions (patches) into columns, and then
   * performs a matrix multiplication with the filters to compute the
   * output maps.
   *
   * Inputs:
   *  - X: Inputs, of shape (N, C*Hin*Win).
   *  - W: Weights, of shape (F, C*Hf*Wf).
   *  - b: Biases, of shape (F, 1).
   *  - C: Number of input channels (dimensionality of input depth).
   *  - Hin: Input height.
   *  - Win: Input width.
   *  - Hf: Filter height.
   *  - Wf: Filter width.
   *  - strideh: Stride over height.
   *  - stridew: Stride over width.
   *  - padh: Padding for top and bottom sides.
   *      For same output height as input, set `padh = (Hf - 1) / 2`,
   *      assuming `strideh = 1`.
   *      More generally, `padh = (Hin*(strideh-1) + Hf - strideh) / 2`
   *      preserves the spatial dimensions of the input.
   *  - padw: Padding for left and right sides.
   *      For same output width as input, set `padw = (Wf - 1) / 2`,
   *      assuming `stridew = 1`.
   *      More generally, `padw = (Win*(stridew-1) + Wf - stridew) / 2`
   *      preserves the spatial dimensions of the input.
   *
   * Outputs:
   *  - out: Outputs, of shape (N, F*Hout*Wout).
   *  - Hout: Output height.
   *  - Wout: Output width.
   */
  N = nrow(X)
  F = nrow(W)
  Hout = as.integer(floor((Hin + 2*padh - Hf)/strideh + 1))
  Wout = as.integer(floor((Win + 2*padw - Wf)/stridew + 1))

  # Create output volume
  out = matrix(0, rows=N, cols=F*Hout*Wout)

  # Convolution - im2col implementation
  parfor (n in 1:N) {  # all examples
    Xn = matrix(X[n,], rows=C, cols=Hin*Win)  # reshape

    # Pad image
    Xn_padded = util::pad_image(Xn, Hin, Win, padh, padw, 0)  # shape (C, (Hin+2*padh)*(Win+2*padw))

    # Extract local image patches into columns with im2col, of shape (C*Hf*Wf, Hout*Wout)
    Xn_padded_cols = util::im2col(Xn_padded, Hin+2*padh, Win+2*padw, Hf, Wf, strideh, stridew)

    # Convolve patches with filters
    outn = W %*% Xn_padded_cols + b  # shape (F, Hout*Wout)
    out[n,] = matrix(outn, rows=1, cols=F*Hout*Wout)  # reshape
  }
}

backward = function(matrix[double] dout, int Hout, int Wout,
                    matrix[double] X, matrix[double] W, matrix[double] b,
                    int C, int Hin, int Win, int Hf, int Wf,
                    int strideh, int stridew, int padh, int padw)
    return (matrix[double] dX, matrix[double] dW, matrix[double] db) {
  /*
   * Computes the backward pass for a 2D spatial convolutional layer
   * with F filters.
   *
   * This implementation uses `im2col` and `col2im` internally.
   *
   * Inputs:
   *  - dout: Gradient wrt `out` from upstream, of
   *      shape (N, F*Hout*Wout).
   *  - Hout: Output height.
   *  - Wout: Output width.
   *  - X: Inputs, of shape (N, C*Hin*Win).
   *  - W: Weights, of shape (F, C*Hf*Wf).
   *  - b: Biases, of shape (F, 1).
   *  - C: Number of input channels (dimensionality of input depth).
   *  - Hin: Input height.
   *  - Win: Input width.
   *  - Hf: Filter height.
   *  - Wf: Filter width.
   *  - strideh: Stride over height.
   *  - stridew: Stride over width.
   *  - padh: Padding for top and bottom sides.
   *  - padw: Padding for left and right sides.
   *
   * Outputs:
   *  - dX: Gradient wrt `X`, of shape (N, C*Hin*Win).
   *  - dW: Gradient wrt `W`, of shape (F, C*Hf*Wf).
   *  - db: Gradient wrt `b`, of shape (F, 1).
   */
  N = nrow(X)
  F = nrow(W)

  # Create output gradient volumes
  dX = matrix(0, rows=N, cols=C*Hin*Win)
  dW = matrix(0, rows=F, cols=C*Hf*Wf)
  db = matrix(0, rows=F, cols=1)

  # Partial derivatives for convolution - im2col implementation
  parfor (n in 1:N) {  # all examples
    doutn = matrix(dout[n,], rows=F, cols=Hout*Wout)

    # Compute dW
    Xn = matrix(X[n,], rows=C, cols=Hin*Win)  # reshape
    Xn_padded = util::pad_image(Xn, Hin, Win, padh, padw, 0)  # shape (C, (Hin+2*padh)*(Win+2*padw))
    Xn_padded_cols = util::im2col(Xn_padded, Hin+2*padh, Win+2*padw, Hf, Wf, strideh, stridew)
    dW += doutn %*% t(Xn_padded_cols)

    # Compute db
    db += rowSums(doutn)
    
    # Compute dX
    dXn_padded_cols = t(W) %*% doutn  # shape (C*Hf*Wf, Hout*Wout)
    dXn_padded = util::col2im(dXn_padded_cols, C, Hin+2*padh, Win+2*padw, Hf, Wf,
                              strideh, stridew, "add")
    dXn = util::unpad_image(dXn_padded, Hin, Win, padh, padw)
    dX[n,] = matrix(dXn, rows=1, cols=C*Hin*Win)  # reshape
  }
}

init = function(int F, int C, int Hf, int Wf)
    return (matrix[double] W, matrix[double] b) {
  /*
   * Initialize the parameters of this layer.
   *
   * Note: This is just a convenience function, and parameters
   * may be initialized manually if needed.
   *
   * We use the heuristic by He et al., which limits the magnification
   * of inputs/gradients during forward/backward passes by scaling
   * unit-Gaussian weights by a factor of sqrt(2/n), under the
   * assumption of relu neurons.
   *  - http://arxiv.org/abs/1502.01852
   *
   * Inputs:
   *  - F: Number of filters.
   *  - C: Number of input channels (dimensionality of depth).
   *  - Hf: Filter height.
   *  - Wf: Filter width.
   *
   * Outputs:
   *  - W: Weights, of shape (F, C*Hf*Wf).
   *  - b: Biases, of shape (F, 1).
   */
  W = rand(rows=F, cols=C*Hf*Wf, pdf="normal") * sqrt(2.0/(C*Hf*Wf))
  b = matrix(0, rows=F, cols=1)
}