At the low level of the implementation of convolutions and other mathemtical operations, the fast code becomes pretty ugly. Taking my convolution code as an example, there are some fast things which make it ugly.

Tensor CnnUtils::convolution(const Tensor& image,Tensor& kernel,const int xStride,const int yStride,bool padding
#if PROFILING
    ,Timer *parentTimer
#endif
){ 
    ...
    if(kernelDimens[1]==3 &&  kernelDimens[2]==3){
        //unrolled 3x3 version 
        ...
    }
    else if(kernelDimens[2]>=8){
        //Do x_i*k_i for 8 in the same row with AVX
        //Saves gathering
        ...
    }
    else{
        const int originalImgYBound = imHeight-yKernelRadius;
        const int originalImgXBound = imWidth-xKernelRadius;
        const int kernelChildSizes0 = kernelChildSizes[0];
        const int kernelChildSizes1 = kernelChildSizes[1];
        ...
        const __m256i paddedImageOffsets = _mm256_setr_epi32(
            0,
            xStride,
            2 * xStride,
            3 * xStride,
            4 * xStride,
            5 * xStride,
            6 * xStride,
            7 * xStride
        );
        for(int l=0;l<paddedImgDimens0;l++){
            int newY = 0;
            int newX = 0;
            //Precomputing multiplications
            int kernelChannel = l*kernelChildSizes0;
            int paddedImageChannel = l*paddedImageChildSizes0-xKernelRadius; //saving the subtractions
            for(int y=yKernelRadius;y<originalImgYBound;y+=yStride){
                int resultRow = newY*resultChildSizes0;
                //Vectorised
                int x=xKernelRadius;
                for(;x+7*xStride<originalImgXBound;x+=8*xStride){
                    int paddedImageChannelShortct = paddedImageChannel + x; 
                    float* __restrict__ resultPtr = resultData+resultRow+newX;
                    //May already have result from another input channel
                    __m256 acc = _mm256_loadu_ps(resultPtr);
                    //Do each individual kernel element across 8 convolutions at once
                    //e.g. do kernel (0,0) multiplied by image (0,0),(0,3),(0,6) ...
                    for(int j=0;j<kernelDimens1;j++){
                        int kernelRow = kernelChannel + j*kernelChildSizes1;
                        int paddedImageRow = paddedImageChannelShortct + (y+j-yKernelRadius)*paddedImageChildSizes1;
                        const float* __restrict__ paddedImageRowBase = &paddedImageData[paddedImageRow];
                        float *kernelRowBase = kernelData+kernelRow;
                        for(int k=0;k<kernelDimens2;k++){
                            const float kernelVal = *(kernelRowBase+k);
                            const __m256i paddedImageIndices = _mm256_add_epi32(paddedImageOffsets,_mm256_set1_epi32(k));
                            const __m256 R = _mm256_i32gather_ps(paddedImageRowBase,paddedImageIndices,4); 
                            const __m256 K = _mm256_set1_ps(kernelVal);  
                            //Add our result
                            acc = _mm256_fmadd_ps(K,R,acc);
                        }
                    }
                    //Save our result
                    _mm256_storeu_ps(resultPtr,acc);
                    newX+=8;
                }
                //Scalar tail
                for(;x<originalImgXBound;x+=xStride){
                    float sum = 0;
                    int paddedImageChannelShortct = paddedImageChannel + x; 
                    for(int j=0;j<kernelDimens1;j++){
                        int kernelRow = kernelChannel + j*kernelChildSizes1;
                        int paddedImageRow = paddedImageChannelShortct + (y+j-yKernelRadius)*paddedImageChildSizes1;
                        float *kernelEndPtr = &kernelData[kernelRow+kernelDimens2];
                        const float* __restrict__ paddedImageDataPtr = &paddedImageData[paddedImageRow];
                        for(float *kernelDataPtr = &kernelData[kernelRow]
                            ;kernelDataPtr<kernelEndPtr;kernelDataPtr++,paddedImageDataPtr++){
                            sum += (*kernelDataPtr) * (*paddedImageDataPtr);
                        }
                    }
                    resultData[resultRow+newX] += sum; 
                    newX++;
                }
                newX=0;
                newY++;
            }
        }
    }
    ...
    return result;
}

Such as:

• Branching for cases which can be optimised and then hardcoding values
• Using intrinsics with ugly names
• Vectorised instructions which result in a scalar tail
• Explicitly loading reused values into variables

I realise that some of these will be done by the compiler (especially in -O3) but making them explicit ensures this and also makes it performant in lower optimisation levels such as debug builds. I also realise that there are more optimisations that could take place (e.g. the 8*xStride) but I like some level of readability. My code is significantly slower than Pytorch (around 3x slower on the same CPU; the CV will not be hearing about this). It seems like some people spend a lot longer writing even uglier code.

Without a machine learning background, I just thought that randomly intialise weights meant a suitably scaled rand() . After struggling with exploding gradients for a while, I found He intialisation which aims to keep the activations having a constant variance of 1. At first I though it was some fancy unnecessary technique (like some parts of ML are) but it drastically improved my accuracy and stability.

Proper profiling tools looked like a bit of a faff and so I made my own. It's a single header file and it outputs a tree of timings. Each timer owns child timers which forms the tree. Timers can be reused to get a mean time for each section of code. In my code, the preprocessor only puts them in the profiling build and they don't particularly overcomplicate the code.

Tensor CnnUtils::convolution(const Tensor& image,Tensor& prePaddedImage,Tensor& kernel,const int xStride,const int yStride
#if PROFILING
    ,Timer *parentTimer
#endif
){
    #if PROFILING
        Timer *prePaddedConvolutionTimer = nullptr;
        Timer *paddingTimer = nullptr;
        if(parentTimer){
            prePaddedConvolutionTimer = parentTimer->addChildTimer("prePaddedConvolution");
            paddingTimer = prePaddedConvolutionTimer->addChildTimer("padding");
        }
    #endif
    // Padding and stuff
    #if PROFILING
        if(parentTimer) paddingTimer->stop();
    #endif 

    Tensor result = convolution(prePaddedImage,kernel,xStride,yStride,false
    #if PROFILING
        ,prePaddedConvolutionTimer
    #endif
    );
    
    #if PROFILING
        if(parentTimer) prePaddedConvolutionTimer->stop();
    #endif 
    return result;
}

With my YouTube eduction, I didn't quite have the correct idea of convolutions. I thought that a kernel was a 2D matrix which you pass over each input channel individually to produce an output channel. This doesn't allow channels to mix like a proper 3D kernel does. However, they are a lot quicker to compute (by a factor of the number of input channels), and seem to do useful processing. When switching my model over to real convolutions, I couldn't beat the accuracy of my fake convolutions for quite a while. For simple tasks, fake convolutions could have uses.