Fast Discounted Cumulative Sums in PyTorch

Overview

TODO: update this README!

Fast Discounted Cumulative Sums in PyTorch

PyPiVersion PythonVersion PyPiDownloads License License: CC BY 4.0 Twitter Follow

This repository implements an efficient parallel algorithm for the computation of discounted cumulative sums and a Python package with differentiable bindings to PyTorch. The discounted cumsum operation is frequently seen in data science domains concerned with time series, including Reinforcement Learning (RL).

The traditional sequential algorithm performs the computation of the output elements in a loop. For an input of size N, it requires O(N) operations and takes O(N) time steps to complete.

The proposed parallel algorithm requires a total of O(N log N) operations, but takes only O(log N) time steps, which is a considerable trade-off in many applications involving large inputs.

Features of the parallel algorithm:

  • Speed logarithmic in the input size
  • Better numerical precision than sequential algorithms

Features of the package:

  • CPU: sequential algorithm in C++
  • GPU: parallel algorithm in CUDA
  • Gradients computation wrt input
  • Both left and right directions of summation supported
  • PyTorch bindings

Usage

Installation

pip install torch-discounted-cumsum

API

  • discounted_cumsum_right: Computes discounted cumulative sums to the right of each position (a standard setting in RL)
  • discounted_cumsum_left: Computes discounted cumulative sums to the left of each position

Example

import torch
from torch_discounted_cumsum import discounted_cumsum_right

N = 8
gamma = 0.99
x = torch.ones(1, N).cuda()
y = discounted_cumsum_right(x, gamma)

print(y)

Output:

tensor([[7.7255, 6.7935, 5.8520, 4.9010, 3.9404, 2.9701, 1.9900, 1.0000]],
       device='cuda:0')

Up to K elements

import torch
from torch_discounted_cumsum import discounted_cumsum_right

N = 8
K = 2
gamma = 0.99
x = torch.ones(1, N).cuda()
y_N = discounted_cumsum_right(x, gamma)
y_K = y_N - (gamma ** K) * torch.cat((y_N[:, K:], torch.zeros(1, K).cuda()), dim=1)

print(y_K)

Output:

tensor([[1.9900, 1.9900, 1.9900, 1.9900, 1.9900, 1.9900, 1.9900, 1.0000]],
       device='cuda:0')

Parallel Algorithm

For the sake of simplicity, the algorithm is explained for N=16. The processing is performed in-place in the input vector in log2 N stages. Each stage updates N / 2 positions in parallel (that is, in a single time step, provided unrestricted parallelism). A stage is characterized by the size of the group of sequential elements being updated, which is computed as 2 ^ (stage - 1). The group stride is always twice larger than the group size. The elements updated during the stage are highlighted with the respective stage color in the figure below. Here input elements are denoted with their position id in hex, and the elements tagged with two symbols indicate the range over which the discounted partial sum is computed upon stage completion.

Each element update includes an in-place addition of a discounted element, which follows the last updated element in the group. The discount factor is computed as gamma raised to the power of the distance between the updated and the discounted elements. In the figure below, this operation is denoted with tilted arrows with a greek gamma tag. After the last stage completes, the output is written in place of the input.

In the CUDA implementation, N / 2 CUDA threads are allocated during each stage to update the respective elements. The strict separation of updates into stages via separate kernel invocations guarantees stage-level synchronization and global consistency of updates.

The gradients wrt input can be obtained from the gradients wrt output by simply taking cumsum operation with the reversed direction of summation.

Numerical Precision

The parallel algorithm produces a more numerically-stable output than the sequential algorithm using the same scalar data type.

The comparison is performed between 3 runs with identical inputs (code). The first run casts inputs to double precision and obtains the output reference using the sequential algorithm. Next, we run both sequential and parallel algorithms with the same inputs cast to single precision and compare the results to the reference. The comparison is performed using the L_inf norm, which is just the maximum of per-element discrepancies.

With 10000-element non-zero-centered input (such as all elements are 1.0), the errors of the algorithms are 2.8e-4 (sequential) and 9.9e-5 (parallel). With zero-centered inputs (such as standard gaussian noise), the errors are 1.8e-5 (sequential) and 1.5e-5 (parallel).

Speed-up

We tested 3 implementations of the algorithm with the same 100000-element input (code):

  1. Sequential in PyTorch on CPU (as in REINFORCE) (Intel Xeon CPU, DGX-1)
  2. Sequential in C++ on CPU (Intel Xeon CPU, DGX-1)
  3. Parallel in CUDA (NVIDIA P-100, DGX-1)

The observed speed-ups are as follows:

  • PyTorch to C++: 387 times
  • PyTorch to CUDA: 36573 times
  • C++ to CUDA: 94 times

Ops-Space-Time Complexity

Assumptions:

  • A fused operation of raising gamma to a power, multiplying the result by x, and adding y is counted as a single fused operation;
  • N is a power of two. When it isn't, the parallel algorithm's complexity is the same as with N equal to the next power of two.

Under these assumptions, the sequential algorithm takes N operations and N time steps to complete. The parallel algorithm takes 0.5 * N * log2 N operations and can be completed in log2 N time steps if the parallelism is unrestricted.

Both algorithms can be performed in-place; hence their space complexity is O(1).

In Other Frameworks

PyTorch

As of the time of writing, PyTorch does not provide discounted cumsum functionality via the API. PyTorch RL code samples (e.g., REINFORCE) suggest computing returns in a loop over reward items. Since most RL algorithms do not require differentiating through returns, many code samples resort to using SciPy function listed below.

TensorFlow

TensorFlow API provides tf.scan API, which can be supplied with an appropriate lambda function implementing the formula above. Under the hood, however, tf.scan implement the traditional sequential algorithm.

SciPy

SciPy provides a scipy.signal.lfilter function for computing IIR filter response using the sequential algorithm, which can be used for the task at hand, as suggested in this StackOverflow response.

Citation

To cite this repository, use the following BibTeX:

@misc{obukhov2021torchdiscountedcumsum,
  author={Anton Obukhov},
  year=2021,
  title={Fast discounted cumulative sums in PyTorch},
  url={https://github.com/toshas/torch-discounted-cumsum}
}
You might also like...
Model summary in PyTorch similar to `model.summary()` in Keras

Keras style model.summary() in PyTorch Keras has a neat API to view the visualization of the model which is very helpful while debugging your network.

torch-optimizer -- collection of optimizers for Pytorch
torch-optimizer -- collection of optimizers for Pytorch

torch-optimizer torch-optimizer -- collection of optimizers for PyTorch compatible with optim module. Simple example import torch_optimizer as optim

A PyTorch implementation of EfficientNet
A PyTorch implementation of EfficientNet

EfficientNet PyTorch Quickstart Install with pip install efficientnet_pytorch and load a pretrained EfficientNet with: from efficientnet_pytorch impor

The easiest way to use deep metric learning in your application. Modular, flexible, and extensible. Written in PyTorch.
The easiest way to use deep metric learning in your application. Modular, flexible, and extensible. Written in PyTorch.

News March 3: v0.9.97 has various bug fixes and improvements: Bug fixes for NTXentLoss Efficiency improvement for AccuracyCalculator, by using torch i

A collection of extensions and data-loaders for few-shot learning & meta-learning in PyTorch

Torchmeta A collection of extensions and data-loaders for few-shot learning & meta-learning in PyTorch. Torchmeta contains popular meta-learning bench

PyTorch Extension Library of Optimized Scatter Operations

PyTorch Scatter Documentation This package consists of a small extension library of highly optimized sparse update (scatter and segment) operations fo

PyTorch Extension Library of Optimized Autograd Sparse Matrix Operations

PyTorch Sparse This package consists of a small extension library of optimized sparse matrix operations with autograd support. This package currently

Reformer, the efficient Transformer, in Pytorch
Reformer, the efficient Transformer, in Pytorch

Reformer, the Efficient Transformer, in Pytorch This is a Pytorch implementation of Reformer https://openreview.net/pdf?id=rkgNKkHtvB It includes LSH

higher is a pytorch library allowing users to obtain higher order gradients over losses spanning training loops rather than individual training steps.
higher is a pytorch library allowing users to obtain higher order gradients over losses spanning training loops rather than individual training steps.

higher is a library providing support for higher-order optimization, e.g. through unrolled first-order optimization loops, of "meta" aspects of these

Comments
  • Sampling code

    Sampling code

    @zhu-han this repo contains the sampling code I mentioned to you. The "iterative" aspect of it is not needed here, we just treat it as a simple way to sample from a distribution.

    Below, is some icefall code called unsupervised.py, that I was going to use to sample CTC transcripts for use in unsupervised training. I believe sampling is more correct than taking the top-one, and will avoid it collapsing to blank.

    # Some utilities for unsupervised training                                                                                                                                                                            
    
    import random
    from typing import Optional, Sequence, Tuple, TypeVar, Union, Dict, List
    
    import math
    import torch
    from k2 import RaggedTensor, Fsa
    from torch import nn
    from torch import Tensor
    import torch_iterative_sampling
    
    Supervisions = Dict[str, torch.Tensor]
    
    
    def sample_ctc_transcripts_ragged(
            ctc_output: Tensor,
            paths_per_sequence: int,
            modified_topo: bool) -> RaggedTensor:
        """                                                                                                                                                                                                               
          ctc_output: a Tensor of shape (N, T, C), i.e. (batch, time, num_symbols),                                                                                                                                       
                     containing normalized log-probs.                                                                                                                                                                     
          paths_per_sequence: The number of separately sampled paths that are requested                                                                                                                                   
                     per sequence                                                                                                                                                                                         
          modified_topo:  True if the system is using the modified CTC topology where two                                                                                                                                 
                    consecutive instances of a nonzero symbol can mean either one or two                                                                                                                                  
                    copies of the original symbol.                                                                                                                                                                        
                                                                                                                                                                                                                          
        Returns a RaggedTensor, on the same device as `ctc_output`, with shape (N *                                                                                                                                       
        paths_per_sequence, None), where 1st index is batch_idx * paths_per_sequence                                                                                                                                      
        + path_idx and 2nd idx is the position in the token sequence.  The returned                                                                                                                                       
        RaggedTensor will have no 0's... those will have been removed, as blanks.                                                                                                                                         
        """
        (N, T, C) = ctc_output.shape
    
        # The 'seq_len' arg below is something specific to the "iterative" part of                                                                                                                                        
        # torch_iterative_sampling, which has to do with "sampling without replacement";                                                                                                                                  
        # here, we don't really want to do "iterative sampling", we just want to                                                                                                                                          
        # sample from the distribution once.                                                                                                                                                                              
    
        probs = ctc_output.exp()
        sampled = torch_iterative_sampling.iterative_sample(probs,
                                                            num_seqs=paths_per_sequence,
                                                            seq_len=1).to(dtype=torch.int32)
        # `sampled` now has shape:                                                                                                                                                                                        
        # (N, T, paths_per_sequence, 1)                                                                                                                                                                                   
        sampled = sampled.squeeze(3).transpose(1, 2)
        # `sampled` now has shape (N, paths_per_sequence, T)                                                                                                                                                              
    
        # identical_mask is of shape (N, paths_per_sequence, T-1), and                                                                                                                                                    
        # contains True at each position if sampled[n,s,t] == sampled[n,s,t+1].                                                                                                                                           
        identical_mask = sampled[:,:,1:] == sampled[:,:,:-1]
    
        if modified_topo:
            # If we are using the modified/simplified CTC topology, it is possible for                                                                                                                                    
            # two consecutive instances of a nonzero symbol to represent either                                                                                                                                           
            # one symbol or two.  We choose either, with probability 0.5.  I think this                                                                                                                                   
            # is correct, perhaps should check though.                                                                                                                                                                    
            identical_mask = identical_mask and (torch.randn(*identical_mask.shape,
                                                             device=identical_mask.device) > 0.5)
        # The following statement replaces repeats of nonzero symbols with 0, so only the                                                                                                                                 
        # final symbol in a chain of identical, consecutive symbols will retain its                                                                                                                                       
        # nonzero value.                                                                                                                                                                                                  
        sampled[:,:,:-1].masked_fill_(identical_mask, 0)
    
        sampled = sampled.reshape(N * paths_per_sequence, T)
    
        # The shape of ragged_sampled would be the same as `sampled`.. it's regular.                                                                                                                                      
        # if you query it, though, it would come up as (N * paths_per_sequence, None).                                                                                                                                    
        ragged_sampled = RaggedTensor(sampled)
    
        # Remove 0's from the ragged tensor, to keep only "real" (non-blank) symbols.                                                                                                                                     
        ragged_sampled = ragged_sampled.remove_values_leq(0)
    
        # note: you can create the CTC graphs with k2.ctc_graph(ragged_sampled, modified={True,False})                                                                                                                    
        # You can turn into a List[List[int]] with ragged_sampled.tolist().                                                                                                                                               
        return ragged_sampled
    
    
    
    def _test_sample_ctc_transcripts_ragged():
        for device in ['cpu', 'cuda']:
            # simple case.. N = 1, T == 2, C == 3                                                                                                                                                                         
            ctc_output = torch.Tensor( [[[ 0., 1., 0. ], [ 1., 0., 0. ] ],
                                        [[ 1., 0., 0. ], [ 0., 0., 1. ] ]]).to(device=device).log()
            r = sample_ctc_transcripts_ragged(ctc_output, paths_per_sequence=1, modified_topo=False)
            print("r = ", r)
            assert r == RaggedTensor('[[1], [2]]', dtype=torch.int32,
                                     device=device)
    
    
    
        for device in ['cpu', 'cuda']:
            # simple case.. N = 1, T == 3, C == 3, with repeats.                                                                                                                                                          
            # We use modified == False, so the repeats should be removed.                                                                                                                                                 
            ctc_output = torch.Tensor( [[[ 0., 1., 0. ], [0., 1., 0.], [ 1., 0., 0. ] ],
                                        [[ 1., 0., 0. ], [0., 0., 1.], [ 0., 0., 1. ] ]]).to(device=device).log()
            r = sample_ctc_transcripts_ragged(ctc_output, paths_per_sequence=1, modified_topo=False)
            print("r = ", r)
            assert r == RaggedTensor('[[1], [2]]', dtype=torch.int32,
                                     device=device)
    
    
    
    if __name__ == "__main__":
        _test_sample_ctc_transcripts_ragged()
    
    
    opened by danpovey 3
Owner
Daniel Povey
Daniel Povey
Fast, general, and tested differentiable structured prediction in PyTorch

Torch-Struct: Structured Prediction Library A library of tested, GPU implementations of core structured prediction algorithms for deep learning applic

HNLP 1.1k Jan 7, 2023
Training RNNs as Fast as CNNs (https://arxiv.org/abs/1709.02755)

News SRU++, a new SRU variant, is released. [tech report] [blog] The experimental code and SRU++ implementation are available on the dev branch which

ASAPP Research 2.1k Jan 1, 2023
An optimizer that trains as fast as Adam and as good as SGD.

AdaBound An optimizer that trains as fast as Adam and as good as SGD, for developing state-of-the-art deep learning models on a wide variety of popula

LoLo 2.9k Dec 27, 2022
Tez is a super-simple and lightweight Trainer for PyTorch. It also comes with many utils that you can use to tackle over 90% of deep learning projects in PyTorch.

Tez: a simple pytorch trainer NOTE: Currently, we are not accepting any pull requests! All PRs will be closed. If you want a feature or something does

abhishek thakur 1.1k Jan 4, 2023
null 270 Dec 24, 2022
A lightweight wrapper for PyTorch that provides a simple declarative API for context switching between devices, distributed modes, mixed-precision, and PyTorch extensions.

A lightweight wrapper for PyTorch that provides a simple declarative API for context switching between devices, distributed modes, mixed-precision, and PyTorch extensions.

Fidelity Investments 56 Sep 13, 2022
A PyTorch repo for data loading and utilities to be shared by the PyTorch domain libraries.

A PyTorch repo for data loading and utilities to be shared by the PyTorch domain libraries.

null 878 Dec 30, 2022
Unofficial PyTorch implementation of DeepMind's Perceiver IO with PyTorch Lightning scripts for distributed training

Unofficial PyTorch implementation of DeepMind's Perceiver IO with PyTorch Lightning scripts for distributed training

Martin Krasser 251 Dec 25, 2022
PyTorch framework A simple and complete framework for PyTorch, providing a variety of data loading and simple task solutions that are easy to extend and migrate

PyTorch framework A simple and complete framework for PyTorch, providing a variety of data loading and simple task solutions that are easy to extend and migrate

Cong Cai 12 Dec 19, 2021
Pretrained ConvNets for pytorch: NASNet, ResNeXt, ResNet, InceptionV4, InceptionResnetV2, Xception, DPN, etc.

Pretrained models for Pytorch (Work in progress) The goal of this repo is: to help to reproduce research papers results (transfer learning setups for

Remi 8.7k Dec 31, 2022