<?xml version="1.0" encoding="UTF-8"?><rss xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:content="http://purl.org/rss/1.0/modules/content/" xmlns:atom="http://www.w3.org/2005/Atom" version="2.0" xmlns:cc="http://cyber.law.harvard.edu/rss/creativeCommonsRssModule.html">
    <channel>
        <title><![CDATA[Deep Dives into Computer Science - Medium]]></title>
        <description><![CDATA[Deep dives into machine learning, databases, distributed systems, and other topics - Medium]]></description>
        <link>https://medium.com/deep-dives-into-computer-science?source=rss----2029d38f3ca2---4</link>
        <image>
            <url>https://cdn-images-1.medium.com/proxy/1*TGH72Nnw24QL3iV9IOm4VA.png</url>
            <title>Deep Dives into Computer Science - Medium</title>
            <link>https://medium.com/deep-dives-into-computer-science?source=rss----2029d38f3ca2---4</link>
        </image>
        <generator>Medium</generator>
        <lastBuildDate>Tue, 19 May 2026 12:22:34 GMT</lastBuildDate>
        <atom:link href="https://medium.com/feed/deep-dives-into-computer-science" rel="self" type="application/rss+xml"/>
        <webMaster><![CDATA[yourfriends@medium.com]]></webMaster>
        <atom:link href="http://medium.superfeedr.com" rel="hub"/>
        <item>
            <title><![CDATA[How to design a high-performance neural network on a GPU]]></title>
            <link>https://medium.com/deep-dives-into-computer-science/how-to-design-a-high-performance-neural-network-on-a-gpu-2f7ada309724?source=rss----2029d38f3ca2---4</link>
            <guid isPermaLink="false">https://medium.com/p/2f7ada309724</guid>
            <category><![CDATA[matrix-multiplication]]></category>
            <category><![CDATA[high-performance]]></category>
            <category><![CDATA[nvidia]]></category>
            <category><![CDATA[gpu]]></category>
            <category><![CDATA[machine-learning]]></category>
            <dc:creator><![CDATA[Kiran Achyutuni]]></dc:creator>
            <pubDate>Sat, 11 Jul 2020 04:50:27 GMT</pubDate>
            <atom:updated>2020-07-11T04:50:27.116Z</atom:updated>
            <content:encoded><![CDATA[<figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*BvuSWZ4dJlpuVE07aEweww.jpeg" /></figure><p>GPUs are essential for machine learning. One could go to AWS or Google Cloud to spin up a cluster of these machines with ease. NVIDIA has industry-leading GPU with Tensor cores with the Volta V100 and Ampere A100 accelerators. Which one to pick for the best performance? What configuration? Which works best for your neural network? There are a number of questions that the machine learning architect must deal with in order to design the fastest neural network at the lowest cost. Besides, just using a machine with GPU and Tensor Cores does not guarantee peak performance. So, as a machine learning architect, how should one approach this problem? <em>For sure, you cannot be hardware agnostic</em>. You need to aware of the capabilities of the hardware in order to get the most performance at the lowest cost.</p><blockquote>As a machine learning architect, how should you design the neural network to maximize the performance on a GPU?</blockquote><p>In this article, we will deep dive to understand the levers that the machine learning architect has to maximize performance. In particular, we will focus on <strong><em>matrix-matrix multiplication</em></strong> in this article because it is the most frequent and heavy-duty (O(n³)) mathematical operation in machine learning.</p><p>Let’s start with a simple fully connected 1 hidden layer neural network:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*-DlqSpSMrLGOVEokuEXwdA.png" /><figcaption>Figure 1: Matrix multiplications at each layer of the neural network with the shapes of the matrix multiply at each step indicated in parenthesis. For example, (B, L1) is the shape of the matrix with B rows and L1 columns. MM1, MM2, … MM5 are various matrix-matrix multiplications.</figcaption></figure><p>From the basic neural network, we can see that at layer L2, we perform 3 matrix-matrix multiplications (1 forward, and 2 backward). At layer L1, we perform 2 matrix-matrix multiplications (1 forward, and 1 backward). In fact, at every layer other than the first layer (L1), we perform 3 matrix multiplications. If there are <em>n </em>layers in a neural net, there are <em>3n-1</em> matrix-matrix multiplications, i.e., it increases linearly as the size of the neural network.</p><p>One quick observation is the case when batch size <em>B=1, </em>i.e., we learn one data point at a time. In this case, matrix-matrix degenerates to matrix-vector multiplication. However, in practice, the batch size is never 1. In Gradient Descent, the entire dataset is considered during each learning step, whereas in Stochastic Gradient Descent, a batch B &gt; 1 (but much less than the entire dataset) is considered at each learning step.</p><p>Going forward in this article, let’s focus on a single matrix-matrix multiplication between 2 matrices A and B of dimensions (M, K) and (K, N) respectively resulting in a matrix C of dimension (M, N).</p><p>The dimensions M, N, K are determined by the architecture of the neural network at each layer. For example, in <a href="https://papers.nips.cc/paper/4824-imagenet-classification-with-deep-convolutional-neural-networks.pdf">AlexNet</a>, the batch size is 128 with a few dense layers of 4096 nodes and an output layer with 1000 nodes. This will result in multiplication of (128,4096) and (4096, 1000) matrices. These are pretty large matrices.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/782/1*90lbbHN-thfWExoTSRDUOw.png" /><figcaption>Figure 2. Source: Reference 1. Tiled matrix multiplication</figcaption></figure><p>What does “large” mean and how are these matrices multiplied? By “large”, we mean any matrix that does not fit in memory. Let’s dive deeper into large matrix multiplication. The matrix multiplication we have learnt in textbooks assumes that the matrix fits in memory. In reality, this may not be the case, especially in machine learning. In addition, finely tuned matrix multiplication algorithms must take into account the memory hierarchy in the computer for optimal performance. The workhorse for multiplying matrices that don’t fit in memory is the <em>tiled/blocked matrix multiplication</em> algorithm. In tiled/block matrix multiplication, the matrices are split into smaller tiles/blocks that will fit into memory, and then compute the portion of the resultant product matrix (See figure 2). Figure 3 shows how the tiled/block matrix-multiplication is recursively applied at every level of the memory hierarchy.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/678/1*515Zr7GgL_31BGEGeWETTw.png" /><figcaption>Figure 3: Tiled/Block matrix-matrix multiplication recursively applied through the complete memory hierarchy of an NVIDIA CPU-GPU system. GEMM stands for General Matrix Multiplication. Figure source: AnandTech.</figcaption></figure><p>We won’t get into the precise tiled matrix multiplication algorithms here, the interested reader is referred to <a href="https://dl.acm.org/doi/10.1145/2063384.2063431">this paper</a>. The library routine in <a href="https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms">BLAS</a> for <strong>Ge</strong>neral <strong>M</strong>atrix <strong>M</strong>ultiplication is called <strong>GEMM. </strong><a href="https://docs.nvidia.com/cuda/nvblas/index.html">NVBLAS</a> is the Nvidia implementation of GEMM that takes advantage of the internal GPU architecture and implements the tiled/block matrix multiplication. PyTorch and TensorFlow link to this library on machines with Nvidia GPU. The library does all the heavy lifting for you. But a poorly designed neural network will almost certainly reduce the performance.</p><p>For us, the immediate goal is to figure out the conditions under which matrix-matrix multiplication executes as fast as possible. This is possible only if the GPU is 100% busy and not idling for data tiles. For this, we need to look at the memory hierarchy and how fast data can potentially move through the memory hierarchy levels.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/860/1*cICXwQUjCy3OtesmDttuoQ.png" /><figcaption>Figure 4: Roofline Model</figcaption></figure><p>Memory hierarchy offers key advantages for improving performance: 1) they hide the latency differences between CPU, GPU, memory components, 2) they take advantage of program locality. Also, a well-designed memory hierarchy offers the highest performance at the lowest cost/byte. To keep the GPUs continuously busy, the data tiles have to rapidly fed to the GPUs. This is determined by the data transfer bandwidth and how fast the GPU is processing the data. This performance metric is captured by the <a href="https://en.wikipedia.org/wiki/Roofline_model"><em>ops: bytes ratio</em></a><em> </em>in the <a href="https://dl.acm.org/doi/10.1145/1498765.1498785">Roofline model </a>(Figure 4).<em> </em>Figure 5 shows how to calculate this from the vendor specs. We see that the ops: bytes ratio is 139 for Volta V100, and 416 for Ampere A100. The larger the ops: bytes ratio, more speedup is possible if the computation was memory or arithmetic bound earlier. In other words, a system with higher ops: bytes ratio is more powerful than a lesser one. This is why Ampere A100 is more powerful than Volta V100.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*8NzkWnx9DbFSAf0-vZZK3g.png" /><figcaption>Figure 5: Calculating the ops: bytes ratio from the vendor specs.</figcaption></figure><p>What does the ops: bytes ratio mean for machine learning and matrix multiply? To see this, we must now look at the compute and data requirements of matrix multiply. <a href="https://en.wikipedia.org/wiki/Roofline_model">Arithmetic intensity</a> is defined as the ratio of the floating-point operations/sec to the bytes. Figure 6 shows how to compute the arithmetic intensity.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*sa1jWzXruMOlfUVVBq7UZg.png" /><figcaption>Figure 6: Computing the arithmetic intensity for matrix multiplication</figcaption></figure><blockquote>If arithmetic intensity &gt; ops:bytes, then matrix multiplication is arithmetic bound else it is memory bound.</blockquote><p>Therefore, the first take away is that the dimensions of the matrices involved should be designed such that the arithmetic intensity is greater than ops: bytes. This will ensure that the GPU is fully utilized. For example, if batch size = 512, N=1024, and M=4096, the arithmetic intensity will be 315, which is greater than 139 for the Volta V100 GPU. Therefore, this matrix multiplication is arithmetic bound on Volta V100 and the GPU will be fully utilized. Figure 7 shows the arithmetic intensities for some common operations in machine learning. The second row corresponds to batch size = 1. The linear layer becomes memory bound and not arithmetic bound in this case. This is the reason why a batch size of 1 is generally not used in production machine learning algorithms.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*ajuyD3oENYbi1sKwn_0l8A.png" /><figcaption>Figure 7. Arithmetic intensities for some common operations in machine learning. Source: NVIDIA documentation. See Reference 1.</figcaption></figure><p>However, ensuring the right arithmetic intensity by choosing the right matrix dimensions is not sufficient to achieve peak arithmetic performance which requires keeping all the tensor cores busy as well. In order to effectively use Nvidia’s Tensor Cores, M, N, K must be a multiple of 8 for FP16 arithmetic or multiple of 16 for FP32 arithmetic. Nvidia core libraries check the dimensions of the matrices and if the conditions are satisfied, then the operations are routed to the Tensor cores. This can result in a 6x speedup on Volta using Tensor cores as compared to using without the Tensor cores. Therefore, the second takeaway is that if the dimensions are not a multiple of 8 or 16, then it is recommended that the dimensions are padded appropriately.</p><p>In your quest to speed up the performance as a machine learning architect, you will invariably face the decision whether to upgrade from Volta to Ampere and pay a higher cost. To this end, you must determine if your neural network is either arithmetic bound or memory bound using the Roofline model. If it is neither, then there is no value in upgrading to a more powerful machine. This is the third takeaway. Nvidia offers tools such as <a href="https://developer.nvidia.com/nsight-compute">Nsight Compute</a> to perform application analysis.</p><p>To summarize:</p><ul><li>matrix-matrix multiply is the most frequent operation in neural network training and inference. The number of matrix multiplications is almost <em>3n</em> the number of layers in a neural network. Hence it is important to compute these as fast as possible.</li><li>The matrices are very large in a neural network. Therefore, we will invariably use a GPU to speed up the matrix multiplications. In order to do, we must understand the ops: bytes ratio of the GPU and design the layers in such a way that the arithmetic intensity is greater than the ops: bytes ratio if possible.</li><li>In order to achieve peak arithmetic performance such that all the Tensor cores are used, the dimensions of the matrices must also meet the requirements imposed by the NVIDIA architecture to use the Tensor core. Usually, it is a multiple of 8 (for FP16 arithmetic) or 16 (for FP32 arithmetic). It is best to see the documentation to ensure the requirement is met.</li><li>You should determine if the application is memory bound or arithmetic bound. If it is neither, then there is no point in upgrading to a more powerful GPU. Otherwise, we can accelerate further by upgrading.</li></ul><p>Being aware of the hardware capabilities and the requirements it imposes for maximizing performance will help in choosing the matrix dimensions and batch size judiciously. This will lead to a design of the neural network such that the training can be completed in the shortest possible time at the lowest cost.</p><p>This article touched on one aspect of accelerating machine learning. There are other dimensions to this problem as well. We will look at some more in the future. Hope you found this article useful. Happy architecting.</p><p><strong>References</strong></p><ol><li><a href="https://docs.nvidia.com/deeplearning/performance/dl-performance-gpu-background/index.html#gpu-arch">Nvidia docs: Deep Learning performance documentation</a></li><li><a href="https://www.anandtech.com/show/12673/titan-v-deep-learning-deep-dive/3">Nvidia: It’s all about Tensor Cores</a>, AnandTech Blog 2018</li><li><a href="https://dl.acm.org/doi/10.1145/2063384.2063431">Fast Implementation of DGEMM on Fermi GPU</a>, Proceeding of the 2011 International Conference on High Performance, Storage, and Analysis</li><li><a href="https://dl.acm.org/doi/10.1145/1498765.1498785">Roofline: An Insightful Visual Performance Model for Floating-Point Programs and Multicore Architectures,</a> <em>Comm. of the ACM</em>, April 2009</li><li><a href="https://developer.download.nvidia.com/video/gputechconf/gtc/2019/presentation/s9624-performance-analysis-of-gpu-accelerated-applications-using-the-roofline-model.pdf">Performance Analysis of GPU-Accelerated Applications using the Roofline Model</a>, NVIDIA, 2019</li><li><a href="https://www.codeproject.com/Articles/1169323/Intel-Advisor-Roofline-Analysis">Intel Advisor Roofline Analysis</a></li><li><a href="https://developer.nvidia.com/nsight-compute">Nvidia Nsight Compute</a></li><li><a href="https://www.nvidia.com/content/dam/en-zz/Solutions/Data-Center/a100/pdf/nvidia-a100-datasheet.pdf">Nvidia Ampere A100 Datasheet</a></li><li><a href="https://images.nvidia.com/content/technologies/volta/pdf/tesla-volta-v100-datasheet-letter-fnl-web.pdf">Nvidia Volta V100 Datasheet</a></li></ol><img src="https://medium.com/_/stat?event=post.clientViewed&referrerSource=full_rss&postId=2f7ada309724" width="1" height="1" alt=""><hr><p><a href="https://medium.com/deep-dives-into-computer-science/how-to-design-a-high-performance-neural-network-on-a-gpu-2f7ada309724">How to design a high-performance neural network on a GPU</a> was originally published in <a href="https://medium.com/deep-dives-into-computer-science">Deep Dives into Computer Science</a> on Medium, where people are continuing the conversation by highlighting and responding to this story.</p>]]></content:encoded>
        </item>
        <item>
            <title><![CDATA[Deep dive into the internals of NumPy's linalg.polyfit()]]></title>
            <link>https://medium.com/deep-dives-into-computer-science/deep-dive-into-the-internals-of-numpys-linalg-polyfit-b60156bc5d1e?source=rss----2029d38f3ca2---4</link>
            <guid isPermaLink="false">https://medium.com/p/b60156bc5d1e</guid>
            <category><![CDATA[numpy]]></category>
            <category><![CDATA[machine-learning]]></category>
            <category><![CDATA[linear-algebra]]></category>
            <category><![CDATA[linear-regression]]></category>
            <category><![CDATA[deep-dives]]></category>
            <dc:creator><![CDATA[Kiran Achyutuni]]></dc:creator>
            <pubDate>Thu, 21 May 2020 10:04:53 GMT</pubDate>
            <atom:updated>2020-05-21T10:04:53.140Z</atom:updated>
            <content:encoded><![CDATA[<figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*BvuSWZ4dJlpuVE07aEweww.jpeg" /><figcaption>Image Source: Unsplash: Willian Justen de Vasconcellos</figcaption></figure><p>Today, a linear regression program in machine learning can be written in just a few lines so easily that we can fail to appreciate the depth of the progress in math and computer systems that enables this simplicity for us. Behind the powerful <a href="https://numpy.org/">numpy</a> and <a href="https://www.scipy.org/">scipy</a> libraries, a huge amount of math and numerical complexity is hidden from us. In this article, we do a deep dive into the <em>internal implementation</em> of linear regression to get a better understanding of what exactly is going on under the covers. In the process, we will discover a variety of elegant linear algebra techniques, the complexity of calling Fortran functions from python, etc.</p><p>In a linear regression problem, we need to solve <em>Ax = b</em>. Here (A, b) is the set of observed data, and we need to find the parameters <em>x</em> that fits this data. There are innumerable mathematical ways by which we can solve the system of linear equations. We are interested in the <em>internal</em> mathematical methods that underly machine learning packages such as numpy/scipy.</p><p>Let’s start with an example of how numpy.linalg.polyfit() is used. (You can download the python code at Github <a href="https://gist.github.com/kirachy/402ed5062cf21c8b94635fa265b02d25">linear regression using numpy polyfit</a>).</p><iframe src="" width="0" height="0" frameborder="0" scrolling="no"><a href="https://medium.com/media/5c03141b7a1bb00d03bb3651944d6fec/href">https://medium.com/media/5c03141b7a1bb00d03bb3651944d6fec/href</a></iframe><p>Here is the visualization of the result. The “*” is the observed data (or generated in our case) and the line is drawn by the parameters that were computed by polyfit. Clearly, polyfit is doing a good job of linear regression.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*SdX28iZVkA_IO1LKHMWIZA.png" /></figure><p>Most articles on the web stop at this point and don’t go any deeper. However, diving deep and trying to understand what is going on will help us appreciate how the solution was computed. So let’s deep dive into np.linalg.polyfit().</p><p>The first observation we make is that np.linalg.polyfit() is <strong>not</strong> fully implemented in python. The heart of its computation happens in a library that was written in highly optimized FORTRAN code! Let’s see how this is so. The call to polyfit() works its way through the following blocks as shown below. I have cut and paste relevant calls from a variety of files so that you can see how it all stitches together:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/960/1*-ItR9MTfrU4xKHb4Acvluw.png" /></figure><p>In block 2, the call to polyfit() will construct a <a href="https://en.wikipedia.org/wiki/Vandermonde_matrix">Vandermonde matrix</a> via a call to numpy.linalg.polyvander(), a special matrix where the columns are in a geometric progression. Assuming the user wants to fit a polynomial of degree 5, then columns up to <em>x</em>⁵ are constructed. Note that the user only needs to supply a 1-D column of <em>x</em> values and polyfit() generates the Vandermonde matrix. By doing so, it is preparing for a polynomial linear regression for this equation:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/398/1*0qcC5UwE3_12WmlZUdH5vw.png" /></figure><p>The original matrix A has transformed as follows:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/426/1*7iiKKPtsDWZwYbsjCC9e9g.png" /></figure><p>Note that this is a rectangular matrix with <em>m</em> &gt; <em>n</em> (deg). Depending on the number of data observations, <em>m</em> can be very large whereas n<em> </em>is usually a small number. Once this matrix is set up, preparations begin to call the matrix implementation in LAPACK (linear algebra library).</p><p>The LAPACK library was written in Fortran 90 way back in the 1980s and 1990s but is updated regularly. The latest open-source LAPACK version is 3.9.0 and was released in 2019. It is also optimized by processor vendors such as Intel, Apple etc. The Fortran library is compiled into a dynamic library based on the architecture of the platform and comes with the operating system. In addition, when you install numpy, 2 shared objects (.so) also comes along with it to make the calls to the C API easy. To view these files, go to the &lt;path to your python installation&gt;/lib/python3.7/site-packages/numpy/linalg/ and you will see the following:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/718/1*uMGAb8kAoQt2RmtkMJizmQ.png" /></figure><p>The two important files here:</p><ol><li><em>_umath_linalg.cpython-37m-darwin.so</em>. : This file converts python data structures into C data structures so that the C interface to the LAPACK can be invoked.</li><li><em>lapack_lite.cpython-37m-darwin.so</em>: This file is the C interface to the LAPACK routines. It calls the LAPACK routines and converts the returns values back to python data structures. The actual LAPACK routines are provided by the OS. For example, on Apple macOS, you can find the platform optimized version at /System/Library/Frameworks/vecLib.framework/libLAPACK.dylib</li></ol><p>Let’s resume how the call flow happens.</p><ol><li>The main entry point to computing the least-squares solution is np.linalg.lstsq(), as shown in block 3.</li><li>In block 4 in the above diagram, np.linalg.lstsq() calls the C interface in <em>_umath_linalg.lstsq_m().</em></li><li>Blocks 5–8 are templates to generate the python-C interface. For example, in block 6, <em>call_@lapack_func@</em> becomes <em>call_sgelsd()</em> defined in block 7. In block 7, <em>LAPACK(@lapack_func@)(…)</em> becomes the call to <em>sgelsd_(…)</em>. This is an extern function that is implemented in the libLAPACK.dylib provided by Apple.</li></ol><p>Now depending on the LAPACK method invoked, there are many ways to compute the linear least squares as shown in block 9.</p><ol><li>QR or LQ factorization</li><li>Complete Orthogonal factorization</li><li>Singular Value Decomposition (SVD) factorization</li><li>Divide and conquer SVD</li></ol><p><em>numpy.linalg.lstsq()</em> has chosen to use the divide-and-conquer SVD methods. How do we know? Look in block 6, and you will find these functions that are specified for template generation. So the chosen functions are SGELSD, CGELSD, ZGELSD, and DGELSD.</p><p>How does the divide and conquer SVD work? The high-level description of the algorithm is described in block 9 is part of code comments, but you can look at all the details in the original <a href="http://www.netlib.org/lapack/lawnspdf/lawn88.pdf">paper</a>. It shows exactly how the least-squares solution is computed. Here is the abstract of the paper:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/980/1*HmWWbUstuDkmKVCQ45pMlA.png" /></figure><p>In order to read the above paper, here are core linear algebra concepts that one needs to know: 1) <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition">Singular Value Decomposition</a> (SVD) 2) <a href="https://en.wikipedia.org/wiki/QR_decomposition">QR Factorization</a> 3) <a href="https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse">Moore-Penrose Pseudoinverse matrix</a>. There are a number of good YouTube videos on linear algebra but the one I would recommend is the <a href="https://www.youtube.com/playlist?list=PL49CF3715CB9EF31D">Gibert Strang Lectures on Linear Algebra (MIT)</a>. His <a href="https://www.amazon.com/Introduction-Linear-Algebra-Gilbert-Strang/dp/0980232775">textbook</a> is equally fabulous. (I am capturing all the essential math concepts required for machine learning <a href="https://medium.com/deep-dives-into-computer-science/essential-math-concepts-to-understand-machine-learning-d852c57a7341">here</a>). From the paper, we see that the current implementation for least squares, i.e, SGELSD, in LAPACK is 9 times faster than its predecessor on bidiagonal matrices! This is the reason why as a first step, SGELSD converts the input matrix to a bidiagonal matrix so that we can benefit from the speedups. Once the SVD is computed, it is trivially easy to compute the least-squares solution. We compute the Moore-Penrose pseudoinverse and then derive the least-squares solution from it. The final step is also detailed in the paper.</p><p>In summary, we have traced the complete path of numpy.linalg.polyfit() all the way down to the core Fortran code. We have discovered that it uses important numerical linear algebra techniques such as singular value decomposition. We also looked at the original paper that came up with the algorithms that numpy.linalg.polyfit() uses and found that the algorithms are significantly faster than the prior implementations. To interface Fortran code seamlessly with Python, careful implementation of the C and Python API interfaces were developed so that the process would not crash. High-quality software engineering practices are required to make this robust. In addition, the OS vendors provided highly tuned libraries to make the LAPACK implementation more efficient for their processor architecture. Thus, advances in math, software engineering, computer architectures, and open source have made numpy.linalg.polyfit() possible and fast.</p><p>Hope you enjoyed the deep dive journey into nump.linalg.polyfit().</p><p>References:</p><ol><li><a href="http://www.netlib.org/lapack/">LAPACK — Linear Algebra Package</a>, <a href="http://www.netlib.org/lapack/lug/">LAPACK User’s guide</a>, <a href="https://www.netlib.org/lapack/lug/node27.html">LAPACK Least squares functions</a></li><li><a href="https://github.com/numpy/numpy">Numpy source code on Github</a></li><li><a href="https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-c/top/lapack-routines.html">Developer Reference for Intel® Math Kernel Library 2020 — C</a> Processor vendors like Intel offer optimized libraries for LAPACK and other mathematical libraries.</li><li><a href="https://software.intel.com/content/www/us/en/develop/tools/math-kernel-library.html">Intel Math Kernel Library</a></li><li><a href="https://developer.apple.com/documentation/accelerate/veclib">vecLib: Apple macOS implementation of LAPACK</a></li><li>You can download the source code of LAPACK <a href="http://www.netlib.org/lapack/#_lapack_version_3_9_0">here</a> to see the Fortran files.</li><li><a href="http://www.netlib.org/lapack/explore-html/d0/db8/group__real_g_esolve_gabc655f9cb0f6cfff81b3cafc03c41dcb.html">SGELSD Fortran source code</a></li><li><a href="http://www.netlib.org/lapack/lawnspdf/lawn88.pdf">Efficient Computation of the Singular Value Decomposition with Applications to Least Squares Problems.</a></li></ol><img src="https://medium.com/_/stat?event=post.clientViewed&referrerSource=full_rss&postId=b60156bc5d1e" width="1" height="1" alt=""><hr><p><a href="https://medium.com/deep-dives-into-computer-science/deep-dive-into-the-internals-of-numpys-linalg-polyfit-b60156bc5d1e">Deep dive into the internals of NumPy&#39;s linalg.polyfit()</a> was originally published in <a href="https://medium.com/deep-dives-into-computer-science">Deep Dives into Computer Science</a> on Medium, where people are continuing the conversation by highlighting and responding to this story.</p>]]></content:encoded>
        </item>
        <item>
            <title><![CDATA[Essential math concepts to understand machine learning]]></title>
            <link>https://medium.com/deep-dives-into-computer-science/essential-math-concepts-to-understand-machine-learning-d852c57a7341?source=rss----2029d38f3ca2---4</link>
            <guid isPermaLink="false">https://medium.com/p/d852c57a7341</guid>
            <category><![CDATA[linear-algebra]]></category>
            <category><![CDATA[calculus]]></category>
            <category><![CDATA[probability]]></category>
            <category><![CDATA[machine-learning]]></category>
            <dc:creator><![CDATA[Kiran Achyutuni]]></dc:creator>
            <pubDate>Sun, 10 May 2020 12:54:24 GMT</pubDate>
            <atom:updated>2020-05-13T18:06:09.780Z</atom:updated>
            <content:encoded><![CDATA[<figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*3I8YZJmkEXt7nHJfQuiAhQ.jpeg" /></figure><p>Machine Learning involves lots of mathematics — linear algebra, calculus, and probability and statistics. They work closely with each other in machine learning algorithms. Linear algebra concepts are used to represent and manipulate data and weights. Calculus is used in the iterative optimization algorithms (e.g., stochastic gradient) to arrive at the optimal weights from loss functions. Probability concepts are used to minimize the uncertainty between the optimal model and the real but unknown data generating process.</p><p>Understanding these math concepts and how they complement each other is fundamental to understanding why and how machine learning algorithms work. Having a firm grasp of these concepts will help you in your machine learning career. In this article, I list some of the important concepts in each of the math areas.</p><p>Important Concepts from Linear Algebra</p><iframe src="" width="0" height="0" frameborder="0" scrolling="no"><a href="https://medium.com/media/e30f058a9a27c70230de7c62bd27a0ea/href">https://medium.com/media/e30f058a9a27c70230de7c62bd27a0ea/href</a></iframe><p>Important Concepts from Probability and Statistics</p><iframe src="" width="0" height="0" frameborder="0" scrolling="no"><a href="https://medium.com/media/d54ca4cb9ae6ad5beb078e514bc8734a/href">https://medium.com/media/d54ca4cb9ae6ad5beb078e514bc8734a/href</a></iframe><p>Important Concepts from Calculus</p><iframe src="" width="0" height="0" frameborder="0" scrolling="no"><a href="https://medium.com/media/3ccf4ae53280e7c6858f1ef1b7236345/href">https://medium.com/media/3ccf4ae53280e7c6858f1ef1b7236345/href</a></iframe><p>Here is a sample of some important research papers in machine learning that have made an impact, and the (pre-requisite) math concepts they use. Familiarize yourself with these math concepts and you will be able to read the paper first-hand and appreciate its contributions.</p><iframe src="https://cdn.embedly.com/widgets/media.html?src=https%3A%2F%2Fairtable.com%2Fembed%2FshrmmiBFKa55xwkan&amp;display_name=Airtable&amp;url=https%3A%2F%2Fairtable.com%2FshrmmiBFKa55xwkan&amp;image=https%3A%2F%2Fstatic.airtable.com%2Fimages%2Foembed%2Fairtable.png&amp;key=a19fcc184b9711e1b4764040d3dc5c07&amp;type=text%2Fhtml&amp;schema=airtable" width="800" height="533" frameborder="0" scrolling="no"><a href="https://medium.com/media/d505c1b36ac20b639a171d0a0849a0d7/href">https://medium.com/media/d505c1b36ac20b639a171d0a0849a0d7/href</a></iframe><p>(<em>Note: I intend to update this article regularly with more information. So please revisit frequently to be up-to-date</em>).</p><p>Here are some excellent books to refresh your math concepts:</p><ol><li><a href="https://mml-book.github.io/book/mml-book.pdf">Mathematics for Machine Learning</a> by Marc Peter Deisenroth, A. Aldo Faisal, Cheng Soon Ong, To be published by Cambridge University Press</li><li><a href="https://www.amazon.com/Introduction-Linear-Algebra-Gilbert-Strang/dp/0980232775">Linear Algebra and its Applications</a>, Gilbert Strang (MIT)</li><li><a href="https://www.amazon.com/Introduction-Probability-2nd-Dimitri-Bertsekas/dp/188652923X">Introduction to Probability</a>, Bertsekas and Tsitsiklis</li></ol><img src="https://medium.com/_/stat?event=post.clientViewed&referrerSource=full_rss&postId=d852c57a7341" width="1" height="1" alt=""><hr><p><a href="https://medium.com/deep-dives-into-computer-science/essential-math-concepts-to-understand-machine-learning-d852c57a7341">Essential math concepts to understand machine learning</a> was originally published in <a href="https://medium.com/deep-dives-into-computer-science">Deep Dives into Computer Science</a> on Medium, where people are continuing the conversation by highlighting and responding to this story.</p>]]></content:encoded>
        </item>
        <item>
            <title><![CDATA[Like backpropagation, is there forward-propagation as well?]]></title>
            <link>https://medium.com/deep-dives-into-computer-science/like-backpropagation-is-there-forward-propagation-as-well-fedb22828b36?source=rss----2029d38f3ca2---4</link>
            <guid isPermaLink="false">https://medium.com/p/fedb22828b36</guid>
            <category><![CDATA[backpropagation]]></category>
            <category><![CDATA[tensorflow]]></category>
            <category><![CDATA[automatic-differentiation]]></category>
            <category><![CDATA[machine-learning]]></category>
            <dc:creator><![CDATA[Kiran Achyutuni]]></dc:creator>
            <pubDate>Mon, 23 Mar 2020 19:04:33 GMT</pubDate>
            <atom:updated>2020-05-13T17:04:26.726Z</atom:updated>
            <content:encoded><![CDATA[<figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*NzyxsrkiLjjyjiIuCf123w.png" /><figcaption>By Berland at English Wikipedia — Transferred from en.wikipedia to Commons by JRGomà., Public Domain, <a href="https://commons.wikimedia.org/w/index.php?curid=4860262">https://commons.wikimedia.org/w/index.php?curid=4860262</a></figcaption></figure><h3>The role of automatic differentiation in machine learning</h3><p>Backpropagation is fundamental to machine learning. It is a technique for computing the gradient of a function. Given there is a technique called <em>back</em>propagation, one wonders if there is a technique called forward-propagation as well? Turns out, there is. Let’s dive into differentiation techniques here to understand which differentiation technique is most relevant for a given machine learning architecture.</p><p>The derivative of a function is a fundamental tool of calculus. There are 3 fundamental ways:</p><ol><li>Symbolic differentiation: Here, we manipulate mathematical functions and produce closed from derivates. The output is what you would see if you were to compute the derivatives by hand. For example, you can visit <a href="https://www.wolframalpha.com/input/?i=derivative+of+x%5E4+sin+x&amp;lk=3">wolfram alpha</a> and see symbolic differentiation in action as follows. The problem with symbolic differentiation is <em>expression swell </em>which can result is high inefficiencies.</li></ol><figure><img alt="" src="https://cdn-images-1.medium.com/max/1024/1*RS2I-sfQ1yugq1Txg0BP4g.png" /></figure><p>2. <a href="https://en.wikipedia.org/wiki/Numerical_differentiation">Numerical differentiation</a>: This method involves the use of finite differences to compute the derivative of a function. This is quite a common method and is used in many areas of engineering such as aerospace engineering. fluid dynamics etc. However, one needs to be careful about round-off errors. Hence, these techniques go to a great extent to avoid these problems. In machine learning, such a high level of precision is not required in computing the gradient since the purpose of machine learning is not optimizing the gradient but optimizing the parameters. So the learning method (such as stochastic gradient descent) can eventually learn the most optimal parameters even if the gradient is not accurate. Also, the complexity of a numerical differentiation method such as <a href="https://en.wikipedia.org/wiki/Euler_method">Euler’s method</a> or <a href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods">Runge-Kutta method</a> is <em>O(n) </em>where<em> n </em>is the number of basic math operations in the function whose derivative needs to be performed (e,g, the function in the neuron such as sigmoid etc). Given that there are a large number of derivatives to be performed in a neural network, an <em>O(n)</em> factor slows down the computation of the gradient. For these reasons, numerical differentiation is not used in machine learning.</p><p>3. <a href="https://en.wikipedia.org/wiki/Automatic_differentiation">Automatic differentiation</a>: In this method, the basic insight is that any <em>arbitrary</em> function coded into a computer program can be broken down to basic mathematical functions such as arithmetic, trigonometric, logarithmic functions. Now any function can be constructed as a computational graph consisting of the basic functions. And the derivative (also a function) can be constructed as a computational graph and the chain rule of calculus can be used to compute the derivative at any node in the graph. Both these graphs can be constructed simultaneously as we build out the neural network. In fact, this is what Tensorflow does. Let’s call these graphs the main graph (MG) and the derivative graph (DG).</p><p>Before we proceed, let’s recap what we are trying to compute. In <a href="https://en.wikipedia.org/wiki/Stochastic_gradient_descent">stochastic gradient descent</a> (SGD), the machine learning algorithm iterates in batches of examples to learn the weights of the network. During each iteration, we need to compute the Jacobian before we proceed to the next iteration. We are trying to compute the Jacobian as follows:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/351/1*_1RtXBQWjFQ9kF7TmVfQbA.png" /></figure><p>Here the neural network has <em>n</em> inputs and <em>m</em> outputs. For example, in linear regression with sigmoid output, <em>m</em> = 1 and <em>n</em> ≥ 1. In image recognition neural net like <a href="https://en.wikipedia.org/wiki/AlexNet">AlexNet</a>, <em>m</em> = 1000 (outputs), and <em>n</em> = 224x224x3 = 150528 (inputs). The gradient in the latter case is a 1000x150528 matrix.</p><p>There are 2 modes for automatic differentiation for computing MG and DG:</p><p>1. Reverse more accumulation: During each iteration of SGD, evaluate MG first. Compute the loss and then compute the DG on the way back to compute the gradient using the loss. Then the weights are adjusted and the next iteration is performed. Computing DG on the way back is exactly the <em>backpropagation</em>.</p><p>2. Forward mode accumulation: Evaluate DG along with MG in each iteration. At the end of the forward step, DG is also available. The gradient can then be applied and the next iteration can be started.</p><p>Do note that with either reverse mode or forward mode, the end result is the same — the value of the gradient computed is the same. But the difference is the cost of computation. Why so? Matrix-matrix operations in the forward accumulation mode are more expensive than matrix-vector operations in the reverse accumulation mode (<em>O(n³)</em> vs <em>O(n²)</em> respectively.).</p><p>As a rule of thumb, if <em>n</em> &gt; <em>m </em>(i.e, inputs are more than outputs),<em> </em>then backpropagation is more efficient. Otherwise, forward mode accumulation would be better. Therefore, while designing your neural net, you must carefully choose your gradient computation algorithm that is most efficient for your application. In many popular neural nets today, the outputs are far less than the inputs, and hence backpropagation is more popular.</p><p>Tensorflow supports both modes of automatic differentiation: <a href="https://www.tensorflow.org/api_docs/python/tf/autodiff/ForwardAccumulator">Tensorflow Forward accumulato</a>r, and <a href="https://www.tensorflow.org/api_docs/python/tf/GradientTape">Tensorflow Gradient Tape (backpropagation</a>). Both are implemented in software. But this is gradually getting pushed to lower levels of computation (see <a href="https://www.intel.ai/fully-automatic-differentiation-for-tensor-expressions/#gs.0t1plz">Intel</a>). There is active research on automatic differentiation (see reference #3 below) and we can expect further advances in the future.</p><p>Hope this article gave you insight into why a specific gradient computation is chosen. For those who want to dive in further, please look at the references below.</p><p>References:</p><ol><li><a href="https://arxiv.org/abs/1502.05767">Automatic Differentiation in Machine Learning: A Survey</a> 2018</li><li><a href="https://arxiv.org/pdf/1811.05031.pdf">A review of Automatic Differentiation and its Efficient Implementation</a> 2019</li><li><a href="https://autodiff-workshop.github.io/2016.html">Autodiff Workshop</a> 2016</li></ol><img src="https://medium.com/_/stat?event=post.clientViewed&referrerSource=full_rss&postId=fedb22828b36" width="1" height="1" alt=""><hr><p><a href="https://medium.com/deep-dives-into-computer-science/like-backpropagation-is-there-forward-propagation-as-well-fedb22828b36">Like backpropagation, is there forward-propagation as well?</a> was originally published in <a href="https://medium.com/deep-dives-into-computer-science">Deep Dives into Computer Science</a> on Medium, where people are continuing the conversation by highlighting and responding to this story.</p>]]></content:encoded>
        </item>
        <item>
            <title><![CDATA[Is there logic behind probability?]]></title>
            <link>https://medium.com/deep-dives-into-computer-science/is-there-logic-behind-probability-e6cca75c3f42?source=rss----2029d38f3ca2---4</link>
            <guid isPermaLink="false">https://medium.com/p/e6cca75c3f42</guid>
            <category><![CDATA[machine-learning]]></category>
            <category><![CDATA[frequentists]]></category>
            <category><![CDATA[logic]]></category>
            <category><![CDATA[probability-theory]]></category>
            <category><![CDATA[bayesian-machine-learning]]></category>
            <dc:creator><![CDATA[Kiran Achyutuni]]></dc:creator>
            <pubDate>Mon, 09 Mar 2020 16:35:38 GMT</pubDate>
            <atom:updated>2020-03-09T16:33:46.435Z</atom:updated>
            <content:encoded><![CDATA[<p>Probability theory plays a fundamental role in machine learning. Why? This is because probability theory concerns reasoning when there is uncertainty. In machine learning, we are given some data to learn from, and we need to reason if the machine learning model that we come up with will function well with <em>uncertain and unseen future</em> data. This seems like a silly question to ask, but can we trust probability theory? Are its foundations solid? Believe it or not, for over 300 years there has been a raging debate on this!</p><p>In 1713, <a href="https://en.wikipedia.org/wiki/Daniel_Bernoulli">Bernoulli</a> first came up with some basic rules to figure out the odds in a game of dice that is rolled repeatedly. This was the start of what is known as “<em>frequentist</em> probability theory”. However, there was a limitation in this — namely, what if you could not repeat the experiment repeatedly? For example, what is the probability of the Earth’s temperature rising by 5 degrees by 2030? We cannot apply the “frequentist probability theory” since we cannot repeat the Earth warming experiment again and again. It’s a one-time event. This limitation was noted by Bayes and Laplace very early on. In 1812, Laplace proposed a set of rules as to what is known as “Bayesian probability theory”. In this approach, one takes into account prior knowledge (Bayes’ law) and can “subjectively” assigns a probability to an event based on the available information at that point in time.</p><p>Assigning probability “subjectively” and taking prior information into account was the bone of contention between the two camps. The frequentists would not accept the “subjectivity” of the Bayesian approach, and the Bayesian’s would find the frequentists theory too limiting to be applied elsewhere. Yet for some strange reason, both the camps were working off of theorems that were eerily identical but derived in different ways, namely the product and sum rules of probability.</p><p>There were 4 important advances in the 20th century that finally put an end to this raging 300-year controversy:</p><ol><li>Kolmogorov, 1933 — <a href="https://www.york.ac.uk/depts/maths/histstat/kolmogorov_foundations.pdf">Foundations of the theory of Probability</a>. He proposed a set of axioms to derive the rules of frequentist probability theory. This is the classical probability theory that we still learn in schools and colleges today. This is easy to understand and hence used to introduce probability to students. Bernoulli’s law is easy to deduce from these axioms. The set of axioms he chose were arbitrary.</li><li>Cox, 1946 — <a href="http://jimbeck.caltech.edu/summerlectures/references/ProbabilityFrequencyReasonableExpectation.pdf">Probability, Frequency, and Reasonable Expectation </a>— published in the American Journal of Physics, 1946. Cox showed that the rules of probability can be derived by applying Boolean logic and calculus. This was the first time it was established that probability can be thought of as an extension of boolean logic. Cox had 2 arbitrary axioms as well.</li><li>Polya, 1954 — <a href="https://www.amazon.com/Mathematics-Plausible-Reasoning-Two-Volumes/dp/1614275572/ref=sr_1_4?keywords=polya&amp;qid=1583751975&amp;sr=8-4">Mathematics and Plausible Reasoning</a>. In this treatise, Polya goes to extensive length in showing how mathematicians and scientists reason when they are in the process of scientific discovery. He came with up with a set of simple rules (<a href="https://en.wikipedia.org/wiki/Syllogism">syllogisms</a>) to explain the thinking process. This treatise has influenced scientists for many generations since.</li><li>Jaynes, 2003 —<a href="https://www.amazon.com/Probability-Theory-Science-T-Jaynes/dp/0521592712/ref=sr_1_1?crid=C3CUUA7FMB9Y&amp;keywords=jaynes+probability+theory+the+logic+of+science&amp;qid=1583752010&amp;sprefix=jaynes+%2Caps%2C382&amp;sr=8-1"> Probability Theory, The Logic of Science</a>. One issue with both Kolmogorov and Cox is that they are based on axioms that are different but arbitrary. So what is it to say that there won’t be a dozen more theories tomorrow for probability theory. So Jaynes does not start with a set of axioms, but instead, he asks what common sense logic do humans follow while reasoning with uncertainty? He calls these the <em>desiderata. </em>These are not axioms but only requirements that any set of axioms must satisfy if it has to make sense for human reasoning. He is inspired by Cox and Polya in this endeavour and then goes about to derive the fundamental rules of probability from the desiderata using Boolean algebra and calculus. And lo and behold, there are only 2 fundamental rules of probability — the product rule and the sum rule. On top of this, he shows that the axioms which Cox and Kolmogorov arbitrarily assumed can also be derived from the product and sum rule! He also shows that there can be no other set of axioms. If they exist, either they can be derived from these two or it must be inconsistent.</li></ol><p>Jaynes and Cox thus showed that probability is an extension of logic since it was derivable from boolean algebra. Also, more importantly, the frequentist and Bayesian probability theories have the same roots (logic). This explains why the product and sum rule are exactly the same in both frequentist and bayesian theories. This diagram shows how everything is interlinked.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/960/1*LR45LcMN0Q__OnQ4KmMibQ.png" /></figure><p>Jaynes and Cox showed that it doesn’t matter which camp you are in, both the camps have the same logical origin. However, the Bayesian way of thinking allows probability rules to be applied to a much larger set of problems where experiments cannot be repeated. Hence it is advantageous to use the Bayesian way of thinking.</p><p>So where exactly are we applying Bayesian thinking in machine learning? Let’s take the simple example of linear regression.</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/96/1*BwQvVBUV1cQh1mydHdSVLQ.png" /></figure><p>where <strong>w</strong><em> </em>is the weight vector, <strong>x</strong> is the input example vector, and <em>y</em> is the computed output. Now, we don’t know the exact data process/function <em>f(x)</em> that generated the given data example <em>(y,x)</em>. In machine learning, we are trying to figure out <em>f(x)</em>. Given the weights <em>w</em>, how certain can we about the computed value <em>y</em> to its true value if we use this data point <em>x</em>? This is where Bayesian way of thinking probability comes in. We can impose any probability distribution we like on this computed value. We normally impose the Gaussian distribution because, among all distributions, this captures most uncertainty (maximum entropy) and least prior assumption about how the data was generated. This is written as</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/182/1*Hnm-BdndfFT_g6ZpT_bBRQ.png" /></figure><p>where <em>N </em>is the Gaussian distribution. Do note that we cannot apply the frequentist view here because we are only computing <em>y</em> once for each example <em>x</em> in the training set for a given <em>w</em>. There is no way we can do repeated experiments here. So the only choice we have here is the Bayesian view — i.e., our belief in the amount of uncertainty of <em>y</em> from its true value, which we are capturing by a Gaussian distribution.</p><p>Now, over the entire training set, we multiply the uncertainty of each example by applying the product rule of probability theory to give us the following:</p><figure><img alt="" src="https://cdn-images-1.medium.com/max/156/1*kjiyoB22pSpKjdp1-nlfGA.png" /></figure><p>Once we formulate this, we can compute the Maximum Likelihood Estimate to find the best <strong>w</strong>. I would refer you to any machine learning textbook to see how this is done to solve linear regression.</p><p>To summarize, we have shown that foundations of probability theory are solidly rooted in logic. Both frequentist and Bayesian approaches have the same logical foundations. However Bayesian approach is more general, powerful, and simpler. The Bayesian view is widely used in machine learning.</p><img src="https://medium.com/_/stat?event=post.clientViewed&referrerSource=full_rss&postId=e6cca75c3f42" width="1" height="1" alt=""><hr><p><a href="https://medium.com/deep-dives-into-computer-science/is-there-logic-behind-probability-e6cca75c3f42">Is there logic behind probability?</a> was originally published in <a href="https://medium.com/deep-dives-into-computer-science">Deep Dives into Computer Science</a> on Medium, where people are continuing the conversation by highlighting and responding to this story.</p>]]></content:encoded>
        </item>
    </channel>
</rss>